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SECTION  I 


INTRODUCTION 

General  purpose  { Gi  )  warheads  applied  against  enemy  installations  can 
produce  devastating  effects.  The  best  results  are  obtained  through  maximum 
reliability  of  conventional  explosive  penetration  and  detonation.  However, 
advances  in  target  hardening  technology  have  created  conditions  where  such 
warheads  may  fail  to  meet  this  goal.  Prior  to  fuze  initiation,  metal  contain¬ 
ment  wall  failure  and  fragmentation  of  the  high  explosive  can  result  from  the 
impact  forces  between  the  warhead  and  a  hardened  target  as  shown  in  Figure 
1,  The  detonation  that  is  expected  to  occur  under  totally  confined  conditions 
may  be  reduced  to  an  unsteady,  rapid  deflagration  with  much  less  damage  to  the 
target.  The  principal  motivation  for  the  work  presented  here  is  the  develop¬ 
ment  of  a  model  which  cat.  be  used  to  predict  whether  detonation  will  occur 
within  the  ruptured  warhead  filled  with  fragmented  high  explosives. 

Understanding  the  fluid  mechanics  of  flame  spreading  through  a  porous 
reactive  solid  medium  has  been  the  subject  of  considerable  research  at  the 
University  of  Illinois  at  Urbana-Champai gn.  Under  the  direction  of  Dr.  Herman 
Krier,  much  work  has  been  done  in  the  analysis  of  these  complex  flows.  A 
strong  background  for  the  research  presented  in  this  report  was  provided  by 
the  earlier  efforts  of:  Van  Tassell  and  Krier  (Reference  1);  Krier  and 
Gokhale  (Reference  2);  Krier,  Rajan,  and  Van  Tassell  (Reference  3);  Dimitstein 
(Reference  4);  Krier,  Dimitstein,  and  Gokhale  (Reference  5);  Krier,  Gokhale, 
and  Hughes  (Reference  6);  Krier  and  Kezerle  (Reference  7);  Butler,  Krier,  and 
lembeck  (Reference  8);  and  Krier,  Dahm,  and  Butler  (Reference  9). 
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Figure  1.  Damage  to  GP  Warhead  from  Impact  with  Hardened  Targ 
[Notice  the  assumed  metal  case  failure  and  internal 
explosive  fracture.] 


The  culmination  of  these  previous  studies  has  been  the  development  of  one 
dimensional,  numerical  models  of  the  transient  events  of  deflagration  to 
detonation  transition.  The  computational  representation  solves  the  inviscid 
conservation  equations  using  finite  difference  techniques.  The  codes  devel¬ 
oped  in  References  7-9  can  accurately  simulate  shock  formation  within  the 
fragmented  propellant  or  explosives,  and  are  used  to  depict  DDT-deflagration 
to  detonation  transition.  The  delay  of  DDT  resulting  from  the  loss  of  mass 
through  the  container  walls  was  modeled  by  a  quasi  one-dimensional  analysis  by 
including  sink  terms  in  the  mass  and  momentum  conservation  equations  (Refer¬ 
ence  9).  Significant  detonation  delays  for  a  variety  of  input  conditions  have 
been  calculated  from  this  quasi  one-dimensional  code  (Reference  8,  9). 

However,  the  loss  of  mass,  momentum,  and  energy  through  an  opening  in  the 
casing  is  a  multi-dimensional  phenomenon.  Even  though  the  reaction  front 
progresses  axially  through  the  damaged  warhead,  the  mass  is  basically  ejected 
in  the  radial  direction.  This  turning  of  the  flow  is  accomplished  by  changes 
in  momentum  and  energy  which  cannot  be  described  accurately  by  a  one-dimen¬ 
sional  representation.  A  multi -dimensional  model  is  therefore  needed  to 
provide  a  much  more  accurate  description  of  the  effects  of  mass  loss  prior  to 
the  prediction  of  any  steady  detonation.  The  research  described  in  the  fol¬ 
lowing  chapters  is  intended  to  provide  the  foundation  for  such  an  analysis. 

Three-dimensional  flame  spreading  analyses  are  now  being  investigated  in 
great  detail  by  Oahm  (Reference  10)  and  have  been  reported  by  Markatos  and 
Kirkcaldy  (Reference  11).  The  Dahm  version  simulates  flow  through  a  frag¬ 
mented  solid  explosive  inside  a  container  by  solving  all  the  inviscid  con¬ 
servation  equations.  These  relations  are  a  coupled  set  of  non-1 inear  hyper¬ 
bolic  partial  differential  equations  which  are  solved  using  finite  difference 


technique.  However,  and  as  expected,  the  numerical  scheme  requires  enormous 
amounts  of  computer  memory  due  to  the  fine  grid  resolution  needed  and  other 
intricacies  of  the  problem.  Consequently,  computation  costs  are  very  high. 

The  Markatos  and  Kirkcaldy  model  attempts  to  simulate  flame  spreading  in  a  bed 
of  solid  propellant  grains  used  as  a  gun  charge.  That  study  utilized  many  of 
the  concepts  presented  in  Reference  3,  and  lacks  the  complexities  associated 
with  radial  venting. 

This  study  presents  the  analysis  for  the  two  dimensional  f'.ow  inside  a 
ruptured  bed  of  granulated  high  explosive.  Although  requiring  the  assumption 
of  a  circumferential  crack,  the  model  does  account  for  radial  components  of 
the  flow.  Therefore,  it  will  be  possible  to  describe  the  turning  transitions 
taking  place  as  mass  vents  through  the  rip  in  the  casing  surrounding  the 
explosive  bed.  True  non-planar  burn  initiation  and  early  flame  spreading 
transients  can  also  be  modeled,  but  most  importantly,  the  two  dimensional 
model  features  will  also  include  viscous  gas  effects.  This  creates  a  more 
realistic  flow  simulation  by  satisfying  the  no-slip  velocity  boundary  con¬ 
dition  at  the  gas-wall  interface,  a  process  not  considered  by  either  Dahm 
(Reference  10)  or  Markatos  and  Kirkcaldy  (Reference  11).  This  greatly  com¬ 
plicates  the  differential  equations  representing  momentum  and  energy  conser¬ 
vation,  because  new  second  derivative  terms  must  be  numerically  differenced. 

The  complete  two-dimensional  model  is  presented  in  the  following  sec¬ 
tion.  Included  are  the  necessary  assumptions,  boundary  conditions,  and  con¬ 
servation  equations  for  the  problem.  The  detailed  derivation  of  these  govern¬ 
ing  relations  is  analyzed  in  Section  III.  Section  IV  deals  with  the  manner  in 
which  these  equations  are  solved,  by  describing  the  numerical  integration 
scheme.  The  computer  code,  preliminary  results,  and  topics  for  possible 
future  work  are  presented  in  Section  V. 
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SECTION  II 


FUNDAMENTAL  CONCEPTS  FOR  THE  TWO-DIMENSIONAL  MODEL 


INTRODUCTION 

The  two-dimensional  model  described  below  is  a  mathematical  representa¬ 
tion  of  the  two-phase  (solid-gas)  fluid  mechanics  within  an  impact-damaged 
warhead.  First,  initial  and  boundary  conditions  are  imposed  on  the  problem. 
Then,  partial  differential  equations  and  algebraic  relations  describing  the 
flow  are  applied  to  all  points  in  the  domain.  Finally,  a  numerical  differ¬ 
encing  method  is  utilized  to  solve  the  system.  This  section  will  present  the 
formulation  of  the  two  dimensional  model,  while  a  derivation  of  the  governing 
conservation  equations  is  examined  in  the  next  chapter. 

ASSUMPTIONS 

Before  attempting  the  representation  of  two  dimensional  flame  spreading, 
it  is  necessary  to  make  several  assumptions  about  the  flow  conditions.  These 
are: 


1.  Warhead  damage  occurs  prior  to  fuze  initiation. 

The  model  presented  here  simulates  the  fluid  mechanics  in  a  ruptured  bon* 
only  and  does  not  include  any  treatment  of  the  time- dependent  case  failure. 

The  impact  with  the  target,  case  fracture,  and  fragmentation  of  the  high 
explosive  tal'.e  place  prior  to  time  t-0. 
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2.  A  circumferential  crack  is  produced  by  the  force  of  the  concussion. 

A  ring-shaped  opening  must  be  considered  in  order  to  yield  a  two-dimen¬ 
sional  representation  of  the  damaged  weapon.  While  this  effectively  slices 
the  casing  into  two  portions,  it  must  be  assumed  that  structural  integrity  is 
still  maintained,  and  the  size  of  the  crack  remains  constant. 

3.  The  explosive  particles  are  treated  as  uniform  spheres. 

The  shock  of  the  warhead's  impact  causes  granulation  of  the  solid  explo¬ 
sive.  The  resulting  fragments  have  a  high  surface-to-vol ume  ratio,  and  can  be 
considered  to  be  millimeter  or  sub-millimeter  spheres. 

4.  The  separated  two-phase  analysis  is  used. 

This  formulation  considers  each  component  as  an  individual  fluid,  with 
the  particles  treated  as  a  pseudo-fluid  (Reference  7).  Each  phase  has  its  own 
set  of  conservation  equations,  with  coupling  terms  to  account  for  interactions 
between  the  solid  and  the  gas. 

5.  The  particles  are  uniformly  distributed  at  time  t=0. 

It  is  assumed  that  the  spherical  particles  are  initially  packed  in  a 
modified  face  centered  cubic  arrangement,  with  porosity  (the  ratio  of  gas 
volume  to  total  volume)  typically  ranging  from  0.2  to  0.3.  However,  porosity 
will  change  with  time,  due  to  combustion  and  rapid  pressure  increases  which 
compress  the  granules. 

6.  The  particles  burn  in  compliance  with  an  assumed  ignition  criteria  and 
combustion  rate  law. 
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If  a  particle  exceeds  a  specified  input  energy  (or  temperature),  it 
ignites  and  starts  producing  gas.  The  burning  rate  of  the  explosive  is 
described  by  r  =  bPn  with  b  and  n  being  material  constants.  When  the  poros¬ 
ity  reaches  approximately  0,99,  the  particles  are  assumed  to  be  burned  out  and 
gas  production  ceases.  The  fuse  initiation  at  time  t=0  is  simulated  by  a  zone 
of  particles  with  temperatures  already  high  enough  to  satisfy  the  ignition 
condition. 

7.  The  gas  phase  obeys  a  non-ideal  equation  of  state. 

The  two  dimensional  model  utilizes  a  co-volume  approach  to  account  for 
molecular  interactions  of  the  gas. 

8.  Interphase  heat  transfer  is  accomplished  through  convection  only. 

Energy  transport  due  to  conduction  and  radiation  effects  is  neglected  in 

this  analysis.  Heat  transfer  takes  place  solely  through  convection.  (See 
Reference  7.) 

9.  The  flow  through  the  damaged  warhead  is  generally  laminar. 

The  major  source  of  turbulence  is  the  two-phase  nature  of  the  flow.  Gas 
moves  around  the  spherical  granules  creating  local  disturbances  and  insta¬ 
bilities  in  its  wake.  However,  the  fluid  properties  are  averaged  over  the 
entire  control  volume,  and  the  turbulent  fluctuations  cancel  out. 

DOMAIN  OF  THE  MODEL 

The  domain  of  the  two  dimensional  model  is  prescribed  by  the  geometry  of 
the  damaged  warhead.  The  warhead  is  presented  as  a  right  circular  cylinder 
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having  solid,  closed  ends.  The  circumferential  wall  is  also  imoermeable, 
except  at  the  opening  where  mass  can  be  ejected.  This  crack  in  the  casing 
must  he  ring-shaped  in  order  to  preserve  the  two-dimensional  nature  of  the 
problem. 

The  domain  is  axi symmetric  as  illustrated  in  Figure  2.  The  radial  coor¬ 
dinate  ranges  from  r=0  at  the  center  line  to  r=R  at  the  circumferential  con¬ 
tainer  wall.  The  axial  coordinate  extends  from  the  near  end,  at  z=0,  to  the 
far  end,  at  z=L.  The  crack  begins  at  z=zc  and  its  length  is  & i 

CONSERVATION  EQUATIONS 

The  fundamental  relations  needed  for  the  analysis  of  fluid  motion  in  the 
ruptured  bomb  are  presented  here  and  are  derived  in  Section  III.  These  equa¬ 
tions  which  govern  the  conservation  of  mass,  momentum,  and  energy  must  be 
satisfied  at  all  points  in  the  domain.  Separated  flow  concepts  require  that  a 
complete  set  of  these  relations  be  written  for  each  phase  with  coupling  terms 
added  to  represent  particle-gas  interaction.  If  viscous  gas  effects  are  also 
included,  this  forms  a  system  of  eight  second  order  nonlinear  partial  dif¬ 
ferential  equations. 

The  relations  listed  below  were  derived  in  cylindrical  coordinates  due  to 
the  axi symmetric  geometry  of  the  problem.  The  necessary  interphase  coupling 
terms  (see  Reference  7)  are  represented  by: 

rc  -  the  gas  production  rate  from  the  burning  particles, 

D  -  the  interph3se  drag, 

and  0  -  the  convective  heat  transfer  between  the  two  constituents. 
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Figure  2.  Schematic  Representation  of  the  Two  Dimensional 
Domain  with  a  Ring-Shaped  Opening 
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Subscripts  g  and  p  indicate  the  gas  and  particle  phases  respectively,  while  r 
and  z  denote  radial  and  axial  directions.  Bulk  densities  are  defined  as: 


•>1  =  *Pg  (1) 

P2  =  J1  -  +  )  pp  (2) 

where  the  porosity,  $,  as  stated  earlier,  represents  the  ratio  of  the  gas 
volume  to  the  total  volume.  The  two-dimensional  conservation  equations  are 
(see  Section  III): 


Conservation  of  Mass  (Gas  Phase) 

3t  ~  7  Jr  (plrug}  "  3?  (plwg)  +  rc 
Conservation  of  Hass  (Particle  Phase) 

3t  ~  7  Jr  (p2rup)  ”  Jz  (p2wp)  "  rc 
Conservation  of  Momentum  (Gas  Phase) 


Radial  Direction 


(piO  =  -  —  —  fp.ru'1)  -  (p.u  w  ) 
3t  g'  r  ar  '•1  g'  az  ^1  g  g' 


32u 


+  wg  l?  af  (rug^  +  +  T  aM  r  ar  (rug) 


(3) 


(4) 


(5) 
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+ 


+  r  u 
c  p 


_a 

ar 


(%♦) 


Axial  Direction 


if  (»lV  -  '  f  aF  -  af  (<=lwg) 


ug‘r  ar  ^  ar  >  __Z  3  az  lr  ar  ^  gJ  az 

3Z 


+  r  w  -  D  -  — —  fP  4} 
c  p  z  az  *•  g*; 


(6) 


Conservation  of  Momentum  (Particle  Phase) 


Radial  Direction 


at  (p2up]  =  *r^(p2rup)  “  n  (p2UpV 


-  r  u  +  D  -  [P  (1  -  *)  ] 
c  p  r  ar  1  p'  T  J 


(7) 


Axial  Direction 


St  (l>2wp)  ■  -  r  jp  (»2rup*p)  -  if  (»?"P 


-  rc»p  *  °t  -  if  IV1  •  ♦» 


(S) 
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Conservation  of  Energy  (Gas  Phase) 


r>  o  o **  *  o  dU„  U „ 

+  u  {2  ff — £)2  +  (~^)2  +  r — ^i2]  -  4  r—9.  +  -2  +  _ a> 

g  1  ar-*  *■  r>  *■  az'  J  3  ^  ar  r  az-* 


aw„  au 


or*  dw  O  -v  1  •»  O  U  ^ 

+  f — 1  ! 1)2+  u  -i  f-I-l  (ru  11  +  u  --4 

1  ar  azJ  g  ar  >-r  ar  *•  g-’J  g  ~2 


[1 J  rru  )  +  -—9.]  +!a_i  (r  ^ja.)  +  w  L"j 
3  ar  lr  ar  ^  Q'  az-1  r  ar  ^  ar-1  a  2 


4  ir  -T1  ’>u  )  ♦  — ^M-q-Du 
3  3z  ‘r  jr  '  gJ  3ZJ 1  v  r  c 


2  2 

u  w 


-  Di«p  +  rc  lEgh™  *  -f  *  -|) 


Conservation  of  Energy  (Particle  Phase) 


3?  tp2Ep)  =  -  r  if  tP2nip(Ep  +  ^)]  -  |f  [p2wp(Ef 

P 


♦  0  *  nr„p  *  D2*p 


2  2 
u  tv 


frChem 

rc^p  2  2> 


In  Equations  (9)  and  (10),  the  total  energies  are  defined  as 


CONSTITUTIVE  EQUATIONS  (Closure) 

The  conservation  of  mass,  momentum,  and  energy  in  the  two  dimensional 
model  is  described  by  eight  nonlinear  partial  differential  equations  con¬ 
taining  eleven  unknown  quantities,  namely:  pg,  pp,  Ug,  Up,  Wg,  Wp,  Pg,  Pp,  Tg, 
Tp,  and  +.  Consequently,  three  additional  relations  are  required  to  provide 
closure  for  the  system: 

1.  An  equation  of  state  for  the  gas  phase 

2.  An  equation  of  state  for  the  particle  phase 

3.  A  procedure  for  determining  porosity  as  a  function  of  local  stress 

conditions. 

In  addition,  several  constitutive  relations  are  necessary  in  order  to 
accurately  represent  the  source/sink  terms.  These  include: 

1.  An  equation  determining  the  gas  production  rate  from  the  burning 
exp! os i ve 

2.  An  equation  for  the  gas  viscosity  coefficient 

3.  A  relation  describing  the  interphase  drag 
.  A  gas-particle  convective  heat  transfer  equation. 
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The  specific  relations  which  were  used  in  this  work  were  taken  from  Reference 
8  and  are  described  in  Appendix  A. 

BOUNDARY  CONDITIONS 

The  conservation  relations  presented  in  Section  II  are  nonlinear  partial 
differential  equations.  In  general,  there  are  a  very  large  number  of  entirely 
different  possible  solutions  for  such  a  system*.  Finding  a  unique  solution 
for  these  equations  requires  specific  knowledge  of  the  physical  characteris¬ 
tics  of  the  problem.  This  information  is  provided  by  the  boundary  conditions. 

For  the  two  dimensional  model  shown  in  Figure  3,  five  conditions  must  be 
specified.  These  are: 

1.  At  the  Center  line  (r-0,  0<^  z  _<L) 

Since  the  damaged  warhead  is  depicted  as  a  right  circular  cylinder  with  a 
ring  shaped  opening,  the  domain  is  symmetric  about  the  center  line.  Conse¬ 
quently,  there  can  be  no  flux  across  this  border,  although  the  flow  can  move 
parallel  to  it.  All  radial  velocity  components  are  set  equal  to  zero,  as  are 
all  gradients  taken  across  this  boundary. 


*  For  example,  the  functions 

u  *  x2  -  y2  u  =  ex  cosy  u  3  Ln  (x2  +  y2) 
are  all  possible  solutions  of  the  second  order  Laplace  equation, 

32u/ax2  +  32u/3y2  =  0 


Case  Opening 
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Figure  3.  Boundary  Conditions  Imposed  on  the  Two  Dimensional  Domai 
(0  <  r  <  R  and  0  <  z  <  L) 


2.  At  the  Near  End  (2*0,  0<_  r  _<R) 

It  is  assumed  that  symmetry  conditions  also  prevail  along  this  bound¬ 
ary.  While  this  would  actually  require  the  inclusion  of  an  additional  crack 
located  outside  the  domain,  these  conditions  will  also  be  used  to  represent  a  * 

solid  wall  where  slip  is  allowed.  Flow  can  move  along  the  boundary  but  cannot 
pass  through  it,  since  axial  velocity  components  and  gradients  are  equal  to 
zero. 

3.  At  the  Far  End  (z=L,  0<_  r^R) 

Here  a  radiative  boundary  condition  is  applied,  allowing  all  disturbances 
to  pass  through  the  wall  as  if  it  did  not  exist.  The  solid  end  is  assumed  to 
be  farther  downstream  outside  the  domain. 


4.  At  the  Circumferential  Wall  (r=R,  0<^ z  <zc  zc+az^  z  <L) 

This  is  a  solid  wall,  requiring  the  radial  velocity  component  to  be  set 
equal  to  zero.  However,  the  no-slip  condition  is  also  applied  at  this  bound¬ 
ary.  Note  that  new  the  axial  velocities  are  forced  to  zero  as  well. 

5.  At  the  Opening  (r=R,  zc  <_  z  _<_  zQ  +  azc) 

It  is  very  difficult  to  define  a  boundary  condition  applicable  at  the 


crack  due  to  the  unknown  status  of  the  flow  beyond  the  opening.  The  radial 
velocity  of  the  gas  at  the  hole  is  determined  by  assuming  the  flow  is 
choked.  This  maximum  value  is  located  at  the  center  of  the  crack,  and  a 
linear  profile  is  assumed,  in  order  to  lower  the  velocity  to  zero  at  the  edges 


of  the  opening.  The  radial  component  of  the  particle  velocity  can  then  be 
found  by  applying  the  interphase  drag  constitutive  relation  as  described  by 
Krier,  Dahm,  an;.  Samuelson  (Reference  12).  The  drag  force  on  the  solids 
produces  a  change  in  momentum  which  is  used  to  predict  a  particle  velocity  at 
the  tear. 

The  specific  details  of  how  these  boundary  conditions  are  included  in  the 
two  dimensional  model  are  explained  in  Chapter  Four  which  deals  with  the 
numerical  methods  employed  to  solve  the  problem. 

INITIAL  CONDITIONS 

Initial  conditions  represent  boundaries  with  respect  to  time  and  are 
applied  at  the  start  of  the  simulation.  Initially,  all  velocities  are  set 
equal  to  zero.  Above-ambient  temperatures  and  pressures  are  determined  from 
polynomial  profiles,  with  several  particles  assumed  to  exceed  the  ignition 
temperature  criteria  and  thus  burning.  Initial  gas  density  values  are 
computed  from  the  equation  of  state,  while  the  solid  phase  density  is  an  input 
parameter.  The  porosity  and  particle  radius  are  also  initially  prescribed  and 
assumed  constant  throughout  the  bed  at  time  t=0. 


SECTION  III 


DERIVATION  OF  THE  TWO  DIMENSIONAL  AXISYMMETRIC  CONSERVATION  EQUATIONS 

The  relations  describing  the  conservation  of  mass,  momentum,  and  energy 
in  the  two-dimensional  model  were  presented  in  the  previous  section  as  Equa¬ 
tions  (3)  -  (10).  The  derivation  of  these  governing  equations  is  discussed 
here,  and  follows  much  of  the  logic  used  by  Krier  and  Kezerle  (Reference  7) 
for  the  one-dimensional  flow  process.  But  additional  complexities  due  to  the 
treatment  of  a  viscous  gas  and  the  cylindrical  geometry  of  the  domain  will  now 
be  incorporated  into  the  analysis.  The  resulting  system  consists  of  eight 
second  order,  coupled,  nonlinear  partial  differential  equations. 

1.  CONSERVATION  OF  MASS 
a.  Gas  Phase  Continuity 


Consider  the  control 
volume  shown,  with  gas 
flowing  through  the  control 
surfaces 
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Since 


Rate  of  Mass 
Accumul atlon 


Rate  of  Rate  of  Rate  of  Mass 

Mass  in  -  Mass  out  +  Generation 


Then, 


3(P  V  ) 

ft9  (°gugAlg)r"  ^pgugAlg^r  +  Ar+  (pgwg^2g)z 


-  fp  w  A,  ’ 
1  9  9  2g- 


+  r  v 

z  +  az  c 


(13) 


where:  pg 

= 

gas  density 

Alg 

= 

cross  sectional  area  of  gas-radial  direction 

% 

= 

cross  sectional  area  of  gas-axial  direction 

V9 

= 

volume  of  gas  within  control  volume 

V 

s 

total  volume  of  control  volume 

U9 

= 

radial  velocity  of  gas 

W9 

= 

axial  velocity  of  gas 

rc 

= 

gas  generation  due  to  combustion  of  solid  particles 

From  the  geometry  of  the  control  volume: 


V  =  r  Ar  A 9  az 

A1  =  raeaz  (14) 

A2  =  rar&e 

V 

Now  with  the  definition  of  porosity  t  =  y*' 
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•w- 


Vg  =  $rArA9AZ 
AJg  =  4>rA0AZ 
A2g  =  4>rArA0 


Substitute  back  for  Ajg,  A2g,  Vg,  and  V  to  get 

9  {p  $rArA6Az) 

at  (pgugKA0Az)r-  (Pgug*rA0A2)r+Ar 


+  ^pgwg*rAr&e^z‘  ^pgwg'(>rArAe^2+Az  +  rcrArA9AZ 


Letting  4>Pg  =  and  dividing  by  rArA0AZ  gives 


=  1  f(plU9r)r~  ^Plugr^r+Ar i  ,  f(plwg)  -  (p1wq)z+az, 

r  1  Ar  J  +  L  t I - J  +  rc 


Now  take  the  limit  as  Ar,  az  >  0  to  get  the  final  form  of  the  gas  continuity 
equation 


ZX  -  1  _L  /  ,  _a  , 

at  -  -  r  aF  (plu9r)  *  3?  (pi"g5  +  rc 


b.  Solid  Phase  Continuity 


From  a  control  volume  analysis  similar  to  that  used  for  the  gas  phase,  one  can 
write: 


3(p  VJ 

3t  ~  ^ppupAlp^r~  ^ppup^lpV+Ar 


+  (p  w  A„_).  -  (< 


(19) 


Here 


pp  =  solid  particle  density 

Alp  =  cross  sectional  area  of  particles-radial  direction 

A2p  *  cross  sectional  area  of  particles-axial  direction 

Vp  ~  volume  of  solid  particles  within  control  volume 
Up  =  radial  velocity  of  particles 
wp  =  axial  velocity  of  particles 


Notice  that  the  rcV  term  is  negative  here,  since  the  mass  of  particles 
decreases  as  gas  is  generated.  Since 

V  v  -  vn  vn 

vp  =  v  -  vg  -§■  =  — s]i  *  i  -  Hb  i  -  ♦  ^20) 

Therefore 

vp  =  (l  -  4>)  rAraeaz 

Alp  =  (1  -  *)  raeaz  (21) 

A2p  =  (l  -  ♦)  rAFA0 
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Now  substitute  for  Ajp,  A£p,  Vp  and  V  to  get 
3[p  (1  -  f)  rArASAz] 

- - it -  tppup{1  '  rA0AZ)r*  [Ppup{l  -  *)  rA9Az]r+Ar 

+  [ppwP(1  ‘  *}  rArA9Iz  "  tppwP(1  "  *]  rAr&eWz 

-  rc  TArAQAZ  (22) 

letting  (1  -  $)Pp  -  ^  ar)d  dividing  by  TAfAQAZ  gives, 

*2.  =  I  >2-pr>r  -  <»2VWr.  t  >2"p>,  -  <P2»p>z+«.  . 

3t  r  l  Ar  i  i  az  1  "  rc 

(23) 

Again  taking  the  limit  as  Ar,  az  *  0,  the  final  form  of  the  solid  phase  con¬ 
tinuity  equation  is: 

=  ‘  r  IF  (p2upr)  ■  af  (p2wp)  -  rc  <24> 

2.  CONSERVATION  OF  MOMENTUM 
a.  Gas  Phase  Momentum 

In  order  to  include  the  gas  viscous  effects,  it  will  first  be  necessary 
to  consider  the  distortion  of  a  moving  fluid  volume  in  three  dimensions.  The 
reader  is  referenced  to  any  of  a  number  of  graduate  level  textbooks  on  fluid 
mechanics  for  detail. 


Extension  Strain 


Radial  Direction 


e  d  t 

rr 


dr  +  -i  (|f )  dt  dr  -  dr 


dr 


=  _i  rlfl'i  ,  IE 


'rr  ar  ’■at-' 


ar 


ar  '■at 


(2! 


Azimuthal  Direction 


E00dt  S 


rde  +  3I  (r  |f)  «*t  d©  -  rde 


*  (r  #}  dt 

r  ae  1  at' 


ee 


I  _i  r  r 
r  ae  tr  at' 


1  rlL  1©  +  r  *2q  1 
r  t30  at  aeat' 


u.  +  1  av 
r  r  ae 


(2< 


Axial  Direction 


E2zdt  = 


dz  +  at  (If)  dt  dz  -  dz 


dz 


az  t3t 


/•  az- 
tat 


'zz 


_  _ a  f az\  _  aw 

az  '■jt'  '  az 


(2; 


Shear  Strain:  r  -  9  Strain 

dt  =  i  [TAN 


i  h  (If)  dtdte 

♦  ™->  - ] 


re 


_  1  r_  3  1  au-i 

2  (r  ar  (?)  +  r  30J 


(28) 


9  -  z  Strain 


+  TAN 


dz 


I 


,  I  ri  +  3V , 
e0Z  2  30  3ZJ 


(29) 


> 

I  z  -  r  Strain 


dt  1  r-r.,,-1  3Z  ^3t )  dtdZ 

erzdt  ’  7  (TAN  - d z - 


1  IF  (If)  dtdr 

+  TAN-1  3r-  ■*?_■ - ] 


dr 


.  r3U  .  3W 

Erz  "  (iz  3?J 


(30) 


Since  it  is  assumed  that  the  gas  phase  behaves  as  a  Newtonian  viscous  fluid, 
the  stress  tensor  can  be  written  as: 
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(31) 


Tij  =  ‘Z  weij  -  Xekkaij 


where  \  =  -  ~  y  by  Stoke's  hypothesis  and  6.,  Is  the  Kronecker  delta. 
Consequently  the  components  of  the  stress  tensor  are: 


rr 


.  -  „  12  in  2  A  3jrul+  1  3V  +  3Wn 
M  >•  3r  3  '•r  ar  r  ae  az' J 


96 


=  -„  [ 2  {“  +  i  B)  .  |  ,i  iiXHl  +  1  31  + 

1  '■r  r  ae'  3  Lp  ar  r  39  az'l 


zz 


=  _  r?  U*  _  1  (J-  ilnyl  +  I  3v  awn 

11  1  az  3  '■r  ar  r  jg  az'J 


T  3  T 

re  1er 


M  [r 

M  1  ar  'T>  r  ae 


X9Z  "  Tze 


=  -  v 


r 3V  +  1  3Wi 
l  az  r  ae-l 


t  =  ,  -  raw  SUi 

zr  rz  "  H  Ur  az' 


(32) 


Now  for  the  two-dimensional  model,  the  azimuthal  gradients  and  velocity  coef¬ 
ficient  are  zero.  Hence  the  stress  components  reduce  to: 
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Rate  of 


Rate  of 


Rate  of 


Sum  of  forces  acting  on 


Momentum  =  Momentum 
Accumulation  In 


-  Momentum  +  the  control  volume  including: 
Out  Shear  Forces 

Gas-Particle  Interaction  Forces 
Pressure  Gradient  Forces 
Forces  associated  with  mass 
addition  or  loss 


S (pq^qU -)  2  2 

-  ft  *  (VsV  '  (V/ig  W  *  (WaVr  <WgV 


Z+AZ 


^TrrAl^r  "  ^TrrAl  V+Ar+  2^TeeA3^sin 


+  ^rzVz  -  (TrzA2W  "  DrV 


+  <PgAlg>r  '  (RgAlg W  +  2<PgA3g>  s1"(f >  +  rcupV  <34> 


Introduced  here  are:  Pg  =  gas  phase  pressure 

Dr  =  radial  component  of  interphase  drag 

A3g  =  cross  sectional  area  of  gas-azimuthal  direction 

It  is  assumed  that  the  pressure  gradient  force  acts  through  the  area  of  the 
gas  only.  It  is  also  envisioned  that  combustion  takes  place  at  the  particle 
velocity,  imparting  momentum  to  the  product  gases.  This  is  accounted  for  by 
the  rcupV  term. 

Like  the  earlier  analysis  of  the  geometry  done  in  Equations  (14)  and  (15) 


A-j  =  Ar&2 


(35) 

(36) 


Now  since  a©  is  small,  then  sin  .  Substituting  this  and  the  rela¬ 

tions  in  Equations  (14),  (15),  (35),  and  (36)  into  the  momentum  Equation  (34) 
yiel ds: 

3(p  4>rArA0Azu  )  2  2 

—  at - =  (PgUg^A6Az)r  -  (pgUgK48Az)r+Ar 

+  {pgugwg+rarA0)z  '  (pgugVrArA0W 

+  (VrrA6AZ>r  "  <VrrA6iz>r+Ar 
+  2  (r^rtz)  +  (trzrArA0)r 

-  (TrzrArA6)r+Ar  -  DrrArA0AZ 

+  (PgtrAOAz)r  -  (Pg*rA0Az)r+Ar  +  2  (Pg*ArAz) 

+  rcuprArA0AZ  (37) 


Rewriting  Equation  (37)  in  terms  of  bulk  density  *  $Pp  and  dividing  by 
volume  rArAOAZ  gives: 


3( 


"lV  .  1  h^W.  ,  .WVq'z- 

3t  r  1  Ar  J  <■  Az  J 
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(38) 


I  r [_ Trr r '  r  ~  ^Trrr^r-t-Ar ,  +  * 96  +  rWz^z  "  ^Trz^z+AZi 
r  *•  Ar  J  r  (  az  ‘ 


1  (pfl*r)r  "  (Pa*r^r+Ar  Po* 

-  D  +  -  —3 — - 9 — ^  +  -S_  +  r  u 

r  r  1  Ar  1  r  c  f 


Now  take  the  limit  as  Ar,  az  ■*•  0  and  simplify  the  pressure  terms: 

a(plU9)  .  .l3(n  Jr)  .-i|0  uw  i  .  r\  .  iM.  +  ll£l  i 

at  >•  ar  ^plug  '  az  'plug  g'  ^r  ar  '  rrr'  r  az  * 


D  -  —  (P  a)  t  r  u 
r  ar  '  g?'  c  p 


(39) 


Substituting  the  results  from  (33)  for  the  stress  terms  and  rearranging  yields 
the  final  form  of  the  radial  direction  gas  momentum  equation: 


at 


I  3 

r  ar 


+ 


1  J.  fi  _1 

1  ar  ‘ r  ar 


(plugr)  '  at  (plugv*g)  +  l,g^aF  tp  aF  (rug)J 
3w 

( ru  )  +  -  9 1 }  +  p  u  -  D  — ( P  * ) 

'  g;  azH  c  p  r  ar  '  g?' 


(40) 


(2)  Gas  Phase  Momentum-Axial  Direction 

Again  consider  a  control  volume  but  now  examine  the  forces  acting  axially 
on  it.  (See  figure  on  following  page.) 


a(p  V  w  ) 

3t  "  ^pgugwgAlgV  "  ^pgugwgAlgV+Ar  + 
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\  \ 


pgWgA2g'z"  pgwgA2g^z+AZ+  ^TrzAl^r 


urzAlV+Ar+  ^zzA2^z"  ^7zzA 2^i 


DzV  +  <pgA2gV  <PgA2gW 


+  rcV 


where:  Dz  =  axial  component  of  interphase  drag 

As  done  previously  for  the  radial  momentum  component,  substitute  geomet 
ric  relations  Equations  (14)  and  (15),  divide  by  the  vol  jme  rArABAz,  and  let 
Pi  a  *pg  to  get: 


3(p,w  )  ,  >,WV  (Plu  w  r), 


(piw«),-  (pt*?) 


1  r ^T_rzr^r~  *Trzr  V+Ara  f*TzzV  ^zz^z+i 

r  1  Ar  i  i  Az 

D  4  [V*-  (V>z+«)  4  r  » 

Z  1  AZ  j  *C  p 


Again  let  Ar,  az  ♦  0 


,(plV  1  a  ,  . 

-at 7?r  ('l“g“jr) 


-  0,  -  T7  <P„4)  ♦  r  W 


Now  substitute  in  the  stress  terms  from  Equation  (33)  and  rearrange  to 
get  the  final  form  of  the  axial  direction  gas  momentum  equation: 
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b.  Solid  Phase  Momentum 

The  derivations  of  the  solid  phase  momentum  equations  are  very  similar  to 
those  done  earlier  for  the  gas  phase,  hut  the  complexities  of  viscous  shear 
are  not  included.  Instead,  interpiiase  drag  forces  become  source  terms  while 
combustion  effects  *re  now  .momentum  sinks.  Finally,  assume  that  solid  phase 
pressure  acts  through  the  particle  area  only. 


(1)  Solid  Phase  Momentum-Radial  Direction 

From  a  control  volume  approach  similar  to  that  used  for  the  gas  phase: 


)  7  7 

3t  ”  ^ppupAlp)r"  (ppupAlp V+&r+  ^ppupwpA2p^i 


^ppupwpA2p^z+az+  ^r^  +  ^pAlpV 
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(45) 


with:  Pp  =  particle  phase  pressure 

4s  in  the  geometry  analysis  done  in  Equations  (20)  and  (21) 
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$)  Araz 


(46) 


A3p  =  - 

Now  substitute  for  Alpf  A2p,  A3p,  Vp,  and  V.  Let  sin  (-^|)  «  ^|  to  obtain: 
3 [p  (l-<t»)  raraeazu  ] 

5t - =  [ppup(1-*)rA0Az]r-  [ppu p ( 1  -♦ ) r aqaz ] P+ Ar 
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+  DrfArA6AZ  +  [Pp(l-4)rA0Az]r-  [Pp( 1-*) rA0Az]r+&r 

+  2  [Pp(l-*)APA2]  ^  -  rcuprAra9AZ  (47) 


Rewrite  Equation  (47)  in  terms  of  bulk  density  p?  =  (l-$)Pp  and  divide  by 
raraeaz: 
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Now  take  the  limit  as  Ar,  az  +  0  and  simplify  the  pressure  term  to  get  the 
final  form  of  the  radial  direction  solid  phase  momentum  equation: 


3  (p2u 


£__  =  I  _JL 

3t  r  3r 


^2Upr)  '  37  (p2UpV  -  rcUp  +  Dr  ‘ 


_3 

dr 


[%(!-♦)] 


(49) 


32 


(2)  Solid  Phase  Momentum-Axial  Direction 

Again  using  the  control  volume  approach,  the  axial  momentum  equation  is 


written  as: 
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Substitute  geometric  relations  Equations  (14)  and  (21),  divide  by  the  volume 
rATAeaz,  and  let  p2  =  (l-$)pp  to  get: 
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(51) 


Take  the  limit  as  Ar,  az  +  0  to  arrive  at  the  final  form  of  the  axial  direc¬ 
tion  solid  phase  momentum  equation: 
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3.  CONSERVATION  OF  ENERGY 
a.  Gas  Phase  Energy 


From  the  first  law  of 
thermodynamics  for  a  viscous 
gas  passing  through  the  control 
surfaces,  we  know  that: 


Rate  of  Energy  Rate  of  Energy 

Accumulation  =  In 


Rate  of  Energy 
Out 


Work  done  by  the 
surroundings  including: 
pressure,  viscous  shear, 
and  interphase  drag 

Thus, 


Heat  addition 
+  from 

combustion  or 
convection 
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where: 


Eg  =  total  gas  energy  CygTg  + 

0  =  rate  of  convective  heat  transfer  from  gas  to  solid 

E^chem  _  gas  chemicai  energy  released  during  combustion 


Now  substitute  in  the  geometry  relations  of  Equations  (14)  and  (15)  and  divide 
by  the  volume.  If  pj  =  pg$,  energy  Equation  (53)  becomes: 
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Combine  energy  and  pressure  terms  and  take  the  limit  as  hr,  hi  *  0  to  get: 
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Substitute  Equation  (33)  for  the  stress  terms  and  take  the  derivatives. 
Rearrange  to  get  the  final  form  of  the  gas  energy  equation: 
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b.  Solid  Phase  Energy 

The  solid  phase  energy  equation  is  derived  using  a  control  volume  analy¬ 
sis  similar  to  the  one  used  for  the  gas  phase.  Once  again,  complications  from 
particle  viscosity  are  neglected  as  they  were  in  solid  phase  momentum  equa¬ 
tions. 
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where:  Ep  =  total  particle  energy  Cv  T  + 

EpC  ^  =  Particle  chemical  energy  released  during  combustion 

Substitute  in  the  geometry  Equations  (20)  and  (21)  divide  by  the  volume, 
and  let  P2  =  pp(l-$)  to  get: 
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(58) 
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Combine  pressure  and  energy  terms  and  take  the  limit  as  at,  az  +  0  to  yield 
the  final  form  of  the  solid  phase  energy  equation 
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SECTION  IV 


NUMERICAL  INTEGRATION  TECHNIQUE 

INTRODUCTION 

In  the  previous  sections,  the  mathematical  foundations  of  the  two-dimen¬ 
sional  model  were  derived  and  discussed  in  great  detail.  The  system  of  par¬ 
tial  differential  equations,  constitutive  relations,  initial  conditions,  and 
boundary  conditions  form  the  basis  of  a  well-posed  problem,  indicating  that  a 
solution  is  possible.  However,  analytical  methods  for  solving  this  set  of 
complex,  coupled,  non-linear  partial  differential  equations  do  not  currently 
exist.  Therefore,  a  computational  numeric  method  must  be  utilized  in  order  to 
determine  the  solution  of  the  system.  The  implementation  of  such  a  finite 
differencing  scheme  is  discussed  here. 

DOMAIN  DISCRETIZATION 

Before  a  numerical  integration  technique  can  be  used  to  solve  a  system  of 
partial  differential  equations,  the  domain  of  the  problem  must  be  re¬ 
defined.  The  conti nous  region  is  replaced  by  a  set  of  control  cells,  which 
represent  locations  where  values  of  dependent  variables  are  stored.  These 
discrete  nodes  make  it  possible  to  approximate  derivatives  which  would  ordi¬ 
narily  be  unknown. 

The  two-dimensional  model  utilizes  a  staggered  grid  to  represent  the 
domain  which  was  described  in  Section  II.  In  this  arrangement,  velocity 
components  are  calculated  at  points  that  lie  on  the  edges  of  the  control  cells 
while  scalar  quantities  are  defined  at  the  centers.  As  shown  In  Figure  4, 


O  Radiol  Velocity  (u)  Node 

□  Axial  Velocity  (w)  Node 

+  Pressure,  Temperature,  Density, 

Porosity  Node 

Figure  4.  Typical  Two  Dimensional  Finite  Difference 
Cell  Utilizing  the  Staggered  Grid 
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radial  velocities  are  determined  at  the  sides  that  are  normal  to  the  r  direc¬ 
tion,  and  axial  velocities  are  calculated  at  faces  perpendicular  to  the  z 
di rection. 

Figure  5  illustrates  how  the  staggered  grid  is  implemented  in  the  two- 
dimensional  model.  The  mass  continuity  equations  (Equations  (3)  and  (4)) 
which  are  used  to  determine  and  p2«  and  the  energy  equations  (Equations  (9) 
and  (10))  which  are  used  to  find  Eg  and  Ep,  are  solved  at  the  scalar  nodes. 

The  radial  momentum  equations  (Equations  (5)  and  (7))  which  solve  for  Ug  and 
Up  are  applied  at  the  displaced  radial  velocity  nodes  while  the  axial  momentum 
equations  (Equations  (6)  and  (8))  that  determine  wg  and  wp  are  employed  at  the 
staggered  axial  velocity  nodes.  The  grid  is  .rranged  so  that  the  boundaries 
of  the  domain  pass  through  sites  containing  normal  velocity  components,  al¬ 
though  the  entire  mesh  extends  beyond  these  borders.  These  exterior  nodes 
will  be  shown  to  be  very  useful  in  dealing  with  the  boundary  conditions  in 
finite  difference  form. 

There  are  three  important  advantages  to  be  gained  by  using  a  staggered 
grid  instead  of  a  conventional  one.  First,  it  provides  second-order  accuracy 

i 

for  the  centered  space  finite  differencing  used  here.  Because  flow  informa¬ 
tion  is  known  at  more  locations,  due  to  the  displaced  velocity  nodes,  deriva¬ 
tive  approximations  can  be  made  over  smaller  distance,  as  explained  by  Dahm, 

Samuelson,  and  Krier  (Reference  12).  The  second  advantage  is  that  the  pres¬ 
sure  difference  between  two  adjacent  nodes  becomes  the  natural  driving  force 
for  the  velocity  component  located  between  these  grid  points.  According  to 
Patankar  (Reference  13)  this  feature  eliminates  the  problems  which  can  arise 
in  the  momentum  equation  due  to  a  wavy  pressure  field.  Finally,  as  will  be 
shown,  the  staggered  grid  makes  the  necessary  boundary  conditions  much  easier 
to  implement  in  finite  difference  form. 
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Figure  5.  The  Staggered  Grid  Applied  to  the  Two  Dimensional  Domain 

[Heavy  border  indicates  the  location  of  the  boundaries  while 
dashed  border  represents  the  opening  in  the  container  wall.] 


THE  FINITE  DIFFERENCE  EQUATIONS 
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Once  the  domain  has  been  discretized,  the  continuous  conservation  equa¬ 
tions  (3)  -  (10)  are  cast  into  finite  difference  form  using  the  “leapfrog1' 
scheme  (Reference  14).  This  centered  time,  centered  space,  explicit  method 
evaluates  unsteady  terms  at  time  levels  t  +  At  and  t  -  At,  while  the  fluxes, 
sources,  and  sinks  are  determined  at  the  current  time  t.  Spatial  derivatives 
are  taken  in  a  similar  manner,  by  evaluating  properties  at  grid  sites  adjacent 
to  the  location  in  question  end  dividing  by  the  distance  between  the  nodes. 

The  finite  difference  forms  of  the  conservation  of  mass,  momentum,  and  energy 
are  listed  in  Appendix  B. 

The  leapfrog  scheme  was  chosen  for  the  two  dimensional  model  based  upon 
the  successful  results  reported  by  other  authors.  Kurihara  (Reference  14) 
showed  that  this  method  could  be  used  to  integrate  the  time  dependent  wave 
equation.  Williams  (Reference  15)  applied  this  scheme  to  the  three  dimen¬ 
sional  conservation  equations  for  an  incompressible,  single-phase,  viscous 
fluid.  In  addition,  Dahm  (Reference  10)  is  currently  obtaining  successful 
results  by  employing  this  technique  to  solve  the  two-phase,  reactive  flow 
equations  in  three  dimensions,  a  problem  similar  to  the  one  being  investigated 
here. 


BOUNDARY  CONDITIONS 

When  partial  differential  equations  are  applied  at  a  boundary  of  the 
domain,  special  conditions  must  be  defined  in  order  to  yield  a  unique  solu¬ 
tion.  This  is  also  true  when  dealing  with  the  finite  difference  form  of  the 
equations.  Hence  the  boundary  conditions  must  now  be  stated  iri  a  finite 
difference  configuration. 
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Recall  from  Section  II  that  there  are  five  boundary  conditions  which  must 
be  specified  for  the  two  dimensional  model.  Each  one  is  considered  sepa¬ 
rately: 

1.  At  the  center  line  (Figure  6) 

The  symmetry  conditions  applied  along  this  boundary  specify  that  all 
radial  velocity  components  and  derivatives  in  the  radial  direction  must  vanish 
at  the  center  line.  Consequently,  the  staggered  grid  is  applied  so  that  the 
u-velocity  nodes  are  positioned  on  the  boundary  with  their  values  initialized 
to  zero.  Since  only  the  radial  momentum  equations  are  evaluated  at  these 
points,  halting  the  integration  prior  to  reaching  the  center  line  leaves  the 
u-velocity  components  unchanged.  Their  values  remain  at  zero  and  thus  satisfy 
the  boundary  condition.  Note  that  this  procedure  eliminates  the  problems  of 
evaluating  the  1/r  terms  as  r  0.  The  singularities  which  normally  occur  at 
the  center  line  do  not  have  to  be  dealt  with,  since  no  computations  are  per¬ 
formed  on  the  nodes  at  this  location.  The  values  at  these  grid  points  are 
preset  to  zero  and  then  left  alone. 

It  is  also  necessary  to  be  able  to  evaluate  the  continuity,  energy,  and 
axial  momentum  equations  along  the  boundary,  although  the  nodes  containing  the 
axial  velocities  and  the  scalar  quantities  do  not  lie  directly  on  the  border 
at  r  =  0.  Consequently,  the  symmetry  conditions  are  applied  by  using  reflec¬ 
tion  points  about  the  centerline.  Recall  from  Paragraph  2  and  Figure  5,  that 
the  mesh  of  grid  points  extends  beyond  the  limits  of  the  domain.  The  values 
stored  in  these  exterior  nodes  are  set  equal  to  their  adjacent  counterpart 
directly  across  the  boundary.  Therefore,  all  derivatives  taken  radially  at  r 
=  0  are  equal  to  zero  and  satisfy  the  boundary  condition. 
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Figure  6.  Symmetry  Boundary  Conditions  Applied  at 
the  Center  Line,  r  =  0 


Figure  7.  Implementation  of  Symmetry  Boundary 

Conditions  at  the  Near  End  Wall,  z  =  0 
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It  should  be  pointed  out  that  the  two  dimensional  model  requires  that  the 
values  stored  in  the  reflected  nodes  consistently  match  those  located  within 
the  domain.  This  is  due  to  the  existence  of  the  second  derivative,  viscous 
force  terms  in  the  gas  phase  momentum  equations.  Dahm  (Reference  10)  explains 
how  the  radial  flux  terms  are  automatically  set  to  zero  when  a  staggered  grid 
is  applied.  However,  because  he  is  solving  inviscid  conservation  equations 
containing  only  fluxes,  sources,  and  sinks,  he  can  store  arbitrary  values  in 
his  external  nodes,  since  eventually  they  get  multiplied  by  radial  velocities 
known  to  be  equal  to  zero.  But  viscous  gas  effects  are  due  to  forces  rather 
than  fluxes,  and  the  second  derivatives  are  not  multiplied  by  these  velocity 
components.  Consequently,  the  symmetry  of  the  reflected  nodes  must  be  main¬ 
tained  in  order  to  provide  an  adequate  boundary  condition  for  these  higher 
order  terms. 

2.  At  the  near  end  wall  (Figure  7) 

Symmetry  conditions  are  also  applied  along  this  end  wall.  Here  the  axial 
velocity  components  are  set  equal  to  zero,  as  are  all  gradients  taken  across 
the  boundary.  The  displaced  w-velocity  nodes  lie  directly  on  the  wall  and  are 
Initialized  to  zero.  These  values  remain  unchanged  since  the  axial  momentum 
equations  are  not  applied  at  these  grid  points.  Reflected  external  nodes  are 
again  implemented  to  provide  symmetry  conditions  for  the  second  derivative 
terms. 

3.  At  the  far  end  (Figure  8) 

In  order  to  allow  flow  disturbances  to  propagate  through  this  boundary,  a 
radiative  end  condition  was  utilized.  This  involved  using  a  Taylor  series 
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RADIATIVE 

BOUNDARY 


Figure  8.  Application  of  the  Radiative  Boundary 
Conditions  to  the  Far  End  Wall,  z  =  L 


‘false'  EXTERIOR  NODES 


Figure  9.  No-Slip  Boundary  Condition  Employed  Along 
with  the  Circumferential  Wall,  r  =  R,  for 
0  <  a  <  z_  and  z_  +  Az.  <  z  <  L 


expansion  with  one-sided  finite  differences  to  determine  the  quantities  to  be 
stored  in  the  external  nodes.  These  predicted  results  are  then  used  in  the 
conservation  equations  to  calculate  values  inside  the  domain. 

For  example,  assume  that  “A"  represents  a  generalized  dependent  variable 
which  is  stored  in  either  a  scalar  or  a  velocity  node.  As  illustrated  in 
Figure  8,  the  nodes  at  (j),  (j-1),  (j-2),  and  ( j -3)  are  located  inside  the 
domain  while  (j+1)  is  an  external  grid  point.  Using  a  Taylor  series  expan¬ 
sion,  the  value  of  A  at  the  outside  node  is: 


A  -  A  +  az  —  I  +  5  A  j 

J+1  J  AZ  3Z  'j  2  322  'j 


(60) 


One-sided  differencing  is  used  to  evaluate  the  derivatives  in  this  approxi¬ 
mation: 


3A  ,  3Aj-  4Aj-l*  Aj-2 
3z  'j  2 AZ 


(61) 


2.  2A 5A .  .+  4A.  A .  , 

l-A  |.  *  -J _ l-.V  _ ill 

3z**  J  [ur 


(62) 


s 

& 


which  when  substituted  into  Equation  (60)  yields  the  predicted  value  to  be 
stored  in  the  exterior  node: 


7A 9A.  .+  3A.  A.  , 
A.  _ i- _ J'1  ,  _ ill 
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(63) 
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4.  At  the  circumferential  wall  (Figure  9) 

Recall  that  this  is  a  solid  wall  where  the  no-slip  condition  is  imposed 
upon  the  flow.  Once  again  the  normal,  radial  velocity  nodes  are  positioned  on 
the  boundary  and  the  values  are  intialized  to  zero.  Integration  of  the  radial 
momentum  equation  is  halted  prior  to  reaching  these  points,  in  order  to  leave 
the  u-velocity  components  unchanged  during  the  numerical  simulation. 

However,  the  no-slip  requirement  can  not  be  imposed  through  the  staggered 
grid,  but  must  be  incorporated  into  the  finite  difference  scheme.  While  the 
boundary  conditions  force  the  axial  velocities  to  be  equal  to  zero  at  the 
wall,  the  grid  points  which  store  these  quantities  are  located  elsewhere,  as 
illustrated  in  Figure  9.  Hence  the  axial  velocity  nodes  can  not  Simply  be 
initialized  to  zero  and  then  left  alone.  Instead,  their  values  must  be  deter¬ 
mined  from  the  discretized  z-momentum  equation,  with  special  considerations 
made  for  the  solid  wall  and  the  no-slip  conditions.  This  is  done  by  using 
"lopsided"  finite  differences  for  approximating  the  radial  derivatives  repre¬ 
senting  viscous  gas  effects. 

For  example,  assume  that  the  axial  velocity  gradient  —  must  be  deter- 

d  i 

mined  at  node  (i)  as  shown  in  Figure  9.  At  any  other  boundary,  the  derivative 
could  be  determined  by  utilizing  the  centered  difference  of: 


3w  ,  wi+l"  wi-l 
3r  h  Ar 


(64) 


This  type  of  approximation  is  possible  because  the  velocities  stored  outside 
the  domain  at  (i+1)  are  known  from  symmetry  or  predicted  by  a  Taylor  series. 
However,  the  values  in  the  exterior  nodes  along  the  solid  wail  are  unknown,  so 


center  differencing  can  not  be  applied.  Instead,  an  imaginary  axial  velocity 
node  is  created  on  the  boundary  at  (i  +  V2),  and  a  lopsided  approximation  is 
used  such  that: 


w  ,  .  3wi-  "i-1 

3r  ’i"  3&r 


(65) 


A  derivation  of  this  formula  can  be  found  in  a  variety  of  textbooks  dealing 
with  numerical  integration  methods.  Now  the  no-slip  condition  can  be  imple¬ 
mented  by  setting  the  value  of  + 1/  equal  to  zero  to  get: 


3w 

ar 


-3V  Vi 
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(66) 


Following  a  similar  line  of  reasoning,  second  derivative  values  are  approxi¬ 
mated  at  the  solid  wall  as: 


-12w^+  4w.j  ^ 
3(ar)Z 


(67) 


Viscous  gas  terms,  consisting  of  gradients  evaluated  in  the  axial  direc¬ 
tion,  require  no  special  treatment  at  the  solid  wall.  Centered  space  differ¬ 
encing  is  used  to  approximate  these  derivatives,  since  all  the  necessary  grid 
points  are  located  inside  the  domain. 

The  radial  flux  terms  in  the  conservation  equations  are  satisfied  at  the 
boundary  by  utilizing  "false"  exterior  nodes,  as  described  by  Dahm  (Reference 
10).  Recall  that  grid  points  located  beyond  the  solid  wall  were  not  required 
for  the  evaluation  of  the  gas  viscosity  terms,  since  lopsided  differencing  was 
employed.  In  addition,  when  a  flux  is  determined  across  the  boundary,  the 
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quantities  stored  in  these  exterior  nodes  are  multiplied  by  radial  velocity 
components,  which  are  known  to  be  equal  to  zero.  Consequently,  arbitrary 
values  can  be  stored  in  these  false  exterior  nodes  in  order  to  facilitate  the 
integration  scheme. 
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5.  At  the  opening  (Figure  10) 

The  flow  conditions  imposed  at  the  rip  in  the  circumferential  wall  are 
based  upon  an  assumed  linear  profile  for  the  radial  component  of  the  gas 
velocity.  These  values  are  stored  in  the  boundary  nodes  at  (i+1)  shown  in 
Figure  10: 
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Here  a^  is  the  gas  speed  of  sound,  which  is  a  function  of  local  temperature 
and  pressure  below  the  hole: 
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This  formula  is  derived  from  the  non-ideal  equation  of  state  in  the  paper  by 
Dahm,  Samuelson,  and  Krier  (Reference  12). 

In  order  to  predict  the  values  of  the  particle  velocity  components  normal 
to  the  opening,  the  interphase  drag  constitutive  relation  was  utilized.  It 
was  assumed  that  this  force  acts  to  pull  the  particles  through  the  hole,  such 
that: 
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NODES  PREDICTED  BY  TAYLOR 

SERIES  AMD  ONE-SIDED  DIFFERENCES 


DETERMINED  FROM  LINEAR  PROFILE 


HOLE 

BOUNDARY 


^P2up^+At  =  2At  °r  ( 2*rp i  +  ^p2up^"At 


(70) 


The  reader  is  again  referred  to  Reference  12  for  a  derivation  of  this  rela¬ 
tion. 

The  remaining  scalar  quantities  and  the  axial  velocity  components  are 
determined  by  using  a  Taylor  series  predictor,  similar  to  the  one  derived  for 
the  end  wall  radiative  boundary  conditions.  If  “A"  is  a  generalized  dependent 
variable,  then 


These  predicted  values  are  stored  in  the  exterior  nodes  at  (i+1)  as  shown  in 
Figure  10.  The  conservation  equations  are  then  applied  to  the  grid  points  at 
(i)  using  the  standard  centered  differencing  techniques. 

STABILITY 

The  stability  of  the  finite  difference  scheme  is  always  a  problem  when 
using  numerical  methods  to  integrate  partial  differencing  equations.  Solu¬ 
tions  which  oscillate  or  grow  without  bound  can  develop  due  to  the  necessary 
mathematical  approximations  required  to  discretize  the  system.  But  steps  can 
be  taken  to  minimize  the  effects  of  these  inherent  instabilities.  The  metho¬ 
dology  used  for  the  two  dimensional  model  is  discussed  below. 

Recall  from  Section  III  that  the  explicit,  centered  time,  centered  space, 
leapfrog  technique  was  used  to  cast  the  conservation  equations  into  finite 
difference  form.  In  any  scheme  such  as  this,  a  bound  exists  on  the  maximum 
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size  of  the  time  step.  If  this  limit  is  exceeded,  oscillations  will  soon 
develop  and  the  predicted  solution  will  become  unstable.  The  conditions  which 
determine  this  maximum  time  step  can  be  found  from  the  Von  Neumann  criteria 
which  is  described  by  Richtmeyer  and  Morton  (Reference  16).  However,  the 
complexity  of  the  governing  system  of  equations  makes  a  complete  analysis  of 
this  kind  extremely  unwieldy  and  difficult.  Consequently,  a  trial  and  error 
method  was  utilized  to  find  a  suitable  time  step,  although  a  Von  Newmann 
stability  condition  should  probably  be  derived  in  the  future. 

Another  source  of  instabilities  is  called  time  splitting.  Although  this 
is  a  direct  consequence  of  the  centered  time  differencing  used  in  the  leapfrog 
scheme,  this  method  was  selected  because  Kurihara  (Reference  14)  showed  that 
it  provided  less  damping  of  the  kinetic  energy.  If  "A"  is  a  generalized 
dependent  variable  and  f(r,z)  represents  the  fluxes,  sources,  sinks,  and 
viscous  gas  terms,  centered  time  differencing  can  be  represented  schematically 
as: 


«t+At  .t-At  , 

— m — s  ,<r-z> 


Rearranging  the  terms  gives: 


At+&t  =  Ztt  ftr.z)1  +  At_At 


(72) 


(73) 


Notice  that  oscillations  can  develop  if  the  value  of  At+At  is  less  than  the 
value  of  A^.  These  instabilities  grow  until  eventually,  two  separate  solu¬ 


tions  form  at  the  even  and  odd  time  steps.  However,  this  problem  can  be 
controlled  through  the  use  of  an  intermediate  smoothing  step.  The  time  filter 
developed  by  Robert  (Reference  17)  was  used  for  the  two-dimensional  model: 


(74) 


A*t+At  =  2At  f(r>2)*t  +  At-At 

A1  =  ♦  e(A*t+At  -  2An  ♦  At_At ) 

where  the  starred  terms  represent  unsmoott.ed  quantities,  and  e  is  the  filter 
parameter.  First,  Equation  (74)  is  used  to  predict  values  at  time  t  +  At,  by 
employing  centered  differencing  and  the  unfiltered  quantities  at  time  t. 

Then,  the  smoothing  step  of  (75)  is  implemented  to  produce  a  filtered  value 
at  t.  Finally,  the  iteration  moves  to  the  next  time  level  and  the  process 
repeats.  This  method  has  been  shown  to  provide  excellent  damping  of  the  time 
splitting  oscillations  of  the  leapfrog  scheme  (Reference  18), 

Instability  due  to  the  finite  differencing  of  non-linear  partial  differ¬ 
ential  equations  appears  through  the  mechanism  of  aliasing.  This  misinterpre¬ 
tation  of  the  high  frequency  waves  is  explained  in  the  paper  by  Krier,  Dahm, 
and  Samuel  son  (Reference  12).  Aliasing  is  controlled  by  evaluating  the  con¬ 
servation  equations  as  derived,  in  their  flux  form  (as  opposed  to  the  convec¬ 
tive  form  with  expanded  derivatives)  (Reference  19). 

The  development  of  shocks  in  the  computational  domain  is  another  source 
of  instability  in  the  numerical  integration.  The  extremely  steep  gradients 
produced  in  the  vicinity  of  the  wave  can  cause  explosive  growth  of  the  solu¬ 
tion.  The  viscous  gas  effects  do  help  to  reduce  the  shock  discontinuities, 
however,  the  second  derivatives  cannot  provide  sufficient  damping  of  these 
strong  pertubations.  Therefore,  artificial  diffusion  terms  are  Introduced 
into  the  conservation  equations  in  a  manner  similar  to  the  one  described  by 
Sod  (Reference  20).  In  order  to  ensure  stability,  it  is  necessary  to  evaluate 
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these  quantities  based  upon  values  calculated  at  the  previous  time  step  (Ref¬ 
erence  21).  The  magnitude  of  the  artificial  diffusion  terms  is  determined  by 
numerical  experimentation. 


SECTION  V 


PRELIMINARY  COMPUTATIONS 


THE  COMPUTER  PROGRAM 

The  fluid  dynamics  concepts  and  numerical  integration  techniques  which 
are  described  in  the  previous  chapters  were  incorporated  into  a  FORTRAN-V 
computer  program  called  DDT20.  This  algorithm  features: 

1.  An  external  data  file 

This  enables  a  user  to  vary  the  input  parameters  without  having  to  re¬ 
compile  the  entire  program,  resulting  in  a  lower  overall  computation  cost. 

2.  A  variable  time  step 

For  each  integration,  the  distance  between  two  adjacent  grid  points  is 
divided  by  the  local  sound  speed  in  the  solid  phase.  That  result  is  then 
multiplied  by  an  input  factor  to  yield  the  size  of  the  next  time  increment. 
This  feature  provides  added  stability  and  produces  better  results  when  steep 
gradients  develop.  The  DDT2D  program  automatically  reduces  the  time  step  so 
that  the  rapid  changes  occurring  near  a  disturbance  can  be  accurately  repre¬ 
sented  without  generating  oscillations  in  the  solution. 

3.  A  choice  of  one  or  two  dimensional  burn  initiation 

Ttie  user  can  seloct  from  two  polynomial  profiles  for  the  initial  temper¬ 
ature  distribution  in  the  solid  phase.  One  option  leads  to  a  one  dimensional 
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initiation,  by  producing  an  area  of  burning  particles  along  the  entire  near 
end  wall.  The  alternative  generates  a  localized  Ignition  zone  at  the  r  *  0,  z 
=  0  corner  of  the  domain,  resulting  In  a  two  dimensional  problem. 

4.  The  ability  to  switch  off  the  viscous  gas  effects 

The  viscous  force  and  dissipation  terms  in  the  gas  momentum  and  energy 
equations  can  be  set  to  2ero  by  simply  changing  the  value  of  an  input  para¬ 
meter.  This  feature  is  particularly  useful  when  examining  the  effects  of  gas 
viscosity  on  the  flow  near  the  opening  in  the  outer  wall.  An  identical  but 
inviscid  case  could  be  used  to  highlight  these  aspects  by  providing  con¬ 
trasting  results. 

These  features  make  D0T20  a  versatile  program  which  should  be  easy  to  use.  A 
listing  of  the  current  version  of  this  code  Is  given  in  Appendix  C. 

The  DDT2D  program  requires  a  great  deal  of  computer  memory  due  to  the 
complex  nature  of  the  fluid  dynamics  problem.  At  each  node,  it  is  necessary 
to  store  the  values  of  eight  primary  variables  at  three  time  levels,  as  well 
as  the  values  of  18  secondary  variables  at  a  single  time  level.  For  the  COC 
CYBER-175  computer  currently  in  use  at  the  University  of  Illinois  at  Urbana- 
Champaign,  the  maximum  core  memory  available  is  131,000  words.  After  allowing 
enough  space  for  the  program  itself,  the  remaining  storage  limits  the  size  of 
the  computational  domain  to  under  3000  nodes.  While  this  should  be  sufficient 
for  most  cases,  additional  memory  can  be  obtained  by  using  the  more  expensive, 
off-line  storage  techniques. 
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THE  BASELINE  CASE 

All  the  tests  described  in  this  chapter  were  conducted  for  a  single  set 
of  fluid  characteristics  and  initial  conditions.  These  input  parameters  for 
the  two-dimensional  baseline  case  are  listed  in  Table  1. 

Recall  from  Section  II  that  the  initial  temperature  and  pressure  condi¬ 
tions  were  prescribed  by  polynomial  profiles  imposed  over  the  domain.  The 
ones  utilized  for  the  two-di irons ional  baseline  case  are  shown  in  Figures  11 
and  12.  Notice  that  at  several  points  located  near  r  *  0  ,  z  =  0  the  temper¬ 
atures  satisfy  the  ignition  criteria  and  the  particles  are  assumed  to  be 
burning.  If  a  one-dironsional  initiation  is  desired,  a  similar  set  of  pro¬ 
files  are  used,  with  variation  occurring  only  in  the  axial  direction. 

The  baseline  case  does  not  include  any  mass  loss  through  the  container 
walls.  This  was  done  in  order  to  simplify  the  problem  during  the  computer 
program  evaluation  stage.  When  the  computer  code  can  accurately  simulate 
flame  spreading  in  a  totally  confined  bed,  the  opening  in  the  circumferential 
wall  can  be  re-introduced  into  the  problem. 
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TABLE  l.  BASELINE  CASE  INPUT  PARAMET1 US 


Parameter  Symbol  Nominal  Value 


Length  of  bed 

L 

3.0  cm 

Radius  of  bed 

R 

0.5  cm 

Particle  radius 

rPo 

100  jjm 

Porosity 

4 

0.30 

Particle  density 

pPo 

1.675  g/cm3 

Ambient  gas  temperature 

Tg 

yO 

300  K 

Ambient  particle  temperature 

300  K 

Ambient  gas  pressure 

\ 

100  kPa 

Ambient  particle  pressure 

100  kPa 

Particle  ignition  energy 

Eign 

4.621  x  109  erg/g 

Gas  viscosity 

“So 

1.8  x  10"4  poise 

Prandtl  Number 

Pr 

0.7 

Gas  constant 

R' 

2967.9  erg/K*gmol 

Gas  phase  specific  heat 

s 

1.51  x  10^  erg/g*K 

Particle  phase  specific  heat 

Cvp 

1.51  x  10^  erg/g*K 

Particle  chemical  energy 

j-cnem 

5.74  x  1030  erg/g*K 

Burning  rate  index 

n 

1.0 

Burning  rate  coefficient 

b 

36.8Ccm/s3/[dyne/cni2]n 

Gas  equation  of  state  coefficient 

bl 

4.0  cm3/g 

Axial  space  increment 

M 

1  mm 

Radial  space  increment 

ar 

1  mm 

Time  filter  coefficient 

c 

0.25 
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TEST  COMPUTATIONS 

Before  any  results  can  be  obtained  from  the  two  dimensional  model,  the 
computer  program  should  be  able  to  successfully  reproduce  the  solutions  of 
other  simpler  cases.  The  failure  to  do  so  may  indicate  the  need  for  some 
numerical  fine  tuning  or  variations  of  input  parameters.  (On  the  other  hand, 
the  inability  to  duplicate  previous  results  may  be  symptomatic  of  a  more 
serious  problem  with  the  code.)  For  these  reasons,  comparisons  with  other 
models  are  essential  for  the  development  of  a  viable  two  dimensional  repre¬ 
sentation  of  the  viscous  gas  flow  in  the  damaged  warhead. 

The  first  case  which  was  examined  was  the  one  dimensional  inviscid,  two- 
phase,  reactive  flow  problem  described  by  Butler,  Lembeck,  and  Krler  (Refer¬ 
ence  8).  In  order  to  model  these  flow  characteristics,  viscous  gas  effects 
were  switched  off  and  initiation  was  fixed  to  take  place  in  one  dimension. 

The  resulting  propagation  of  the  flame  front  is  plotted  in  Figure  13.  Notice 
the  excellent  correlation  between  the  three  different  versions  of  the  same 
problem.  The  Dahm  (Reference  10)  three  dimensional  code  (run  with  a  one 
dimensional  initiation)  predicts  almost  identical  flame  front  locations  while 
the  purely  one  dimensional  model  produces  slightly  lower  values.  These 
results  indicate  that  the  DDT2D  computer  code  can  accurately  model  one  dimen¬ 
sional  ,  inviscid  flows. 

However,  this  is  not  the  case  for  the  more  complicated  two  dimensional 
problems.  This  time  the  D0T2D  program  produced  the  locus  of  ignition  fronts 
shown  in  Figure  14.  Note  that  Initially  the  flame  front  location  varies  in 
both  the  axial  and  radial  directions.  As  time  progresses,  the  variations  in 
the  r-direction  become  less  and  less,  until  the  flow  is  nearly  one  dimen¬ 
sional,  after  t  =  1,97  u$,  While  this  result  is  qualitatively  correct,  the 


Figure  14.  Flame  Front  Locus  Showing  the  Transition  from  a  Two-Dimensional 
Initiation  to  One  Dimensional  Flow 


event  takes  place  much  too  fast.  Based  upon  results  obtained  by  Dahm  (Refer¬ 
ence  10),  the  transition  to  one  dimensional  flow  should  not  occur  until  nearly 
50  us.  Obviously,  there  is  too  much  energy  in  the  system  causing  the  parti - 
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cles  to  ignite  prematurely. 

Therefore,  it  was  decided  to  attempt  to  run  a  case  with  a  two  dimensional 
initiation  and  gas  viscosity  effects.  It  was  assumed  that  the  inclusion  of 
the  viscous  dissipation  terms  in  the  gas  energy  equation  would  slow  the 
advance  of  the  flame  front.  Unfortunately,  no  meaningful  results  could  be 
obtained  since  almost  immedi ately,  spatial  oscillations  developed.  These  wild 
swings  in  the  radial  velocity  components  (shown  in  Figure  15)  eventually 
forced  the  program  to  stop  after  11.33  ys.  Respite  numerous  attempts  at 
damping  these  numerical  oscillations  with  adjustments  to  the  input  parameters, 
no  significant  progress  was  made. 
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POSSIBLE  FOLLOW  UP 

Although  the  current  algorithm  used  In  the  DDT?D  code  may  not  be  entirely 
free  from  all  errors,  the  foundations  are  correct  and  work  on  this  two-dimen¬ 
sional,  viscous  gas  model  should  be  continued.  Reasonable  results  have  been 
obtained  for  an  invlscid  one  dimensional  analysis.  For  this  reason,  it  is 
felt  that  this  program  can  eventually  be  used  to  solve  the  two-dimensional 
viscous  flow  problem  without  major  alterations.  Suggested  areas  for  possible 
changes  in  DDT2D  are, 

1.  Provide  better  initial  conditions 

The  initial  temperature  and  pressure  profiles  which  are  currently  used  in 
the  two  dimensional  model  actually  do  not  satisfy  the  conservation  equations 
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at  time  t  =  0.  Consequently,  the  flow  must  quickly  adjust  in  order  to  meet 
the  requirements  of  these  governing  relations.  This  momentary  instability  may 
be  the  trigger  which  touches  off  the  oscillations  which  develop  downstream. 
Therefore,  better  results  may  be  obtained  from  initial  conditions  based  upon 
established  flows  from  inviscid  models,  such  as  a  two-dimensional  version  of 
the  Dahm  code,  given  in  Reference  10. 

2.  Determine  an  optimal  time  step 

It  may  be  advantageous  to  use  a  Von  Neumann  stability  analysis  (Reference 
16)  to  determine  the  formula  for  the  optimal  size  time  step.  This  could  then 
be  incorporated  into  the  computer  program,  enabling  it  to  determine  the  ideal 
time  increment  for  each  integration,  and  insuring  the  stability  of  the 
differencing  scheme. 

3.  Change  the  finite  differencing  scheme 

The  inviscid  conservation  relations  form  a  set  of  hyperbolic  partial 
differential  equations  which  can  be  solved  numerically  by  using  the  leapfrog 
scheme.  However,  when  the  second  derivative  terms  representing  viscous  forces 
are  included,  the  nature  of  the  governing  equations  may  become  parabolic  or 
elliptic.  If  this  is  the  case,  then  a  different  finite  differencing  scheme 
must  be  implemented. 
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APPENDIX  A 


AUXILIARY  DATA  BASE  AND  CONSTITUTIVE  RELATIONS 

The  eight  conservation  equations  presented  in  Section  II  contain  11 
unknown  quantities.  Determining  the  values  of  these  parameters  requires  the 
addition  of  three  state  relations,  thereby  providing  closure  to  the  system. 
Moreover,  constitutive  equations  are  needed  to  describe  the  gas  viscosity  and 
the  interactions  between  phases.  This  appendix  details  these  necessary  sup¬ 
plemental  relations  which  were  used  for  this  report. 

1.  Gas  Phase  Equation  of  State 

A  non-ideal  equation  of  state  was  utilized: 

pg  "  RVg  (1  *  Vg>  (#-'l 

The  value  of  the  coefficient  b^  is  based  upon  a  fit  of  CJ  detonation  data 
produced  by  the  TIGER  code  {Reference  22)  and  is  assumed  to  remain  constant. 

2,  Solid  Phase  Equation  of  State 

The  density  Gf  the  explosive  particles  was  specified  by  applying  the  Tait 
equation  (see  Coyne,  Butler,  and  Krier,  Reference  23): 
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3.  Porosity-Stress  Relation 


There  are  actually  no  suitable  correlations  for  the  dynamic  porosity- 
stress  relation  for  a  packed  bed  of  compressible  solid  particles  over  tha 
range  of  conditions  imposed  by  this  problem.  Consequently,  an  approximate 
numerical  predictor  method  was  used.  At  each  time  step,  an  equilibrium  con¬ 
dition  is  applied,  setting  the  particle  pressure  equal  to  that  of  the  gas. 
Porosity  is  then  determined  by  dividing  the  particle  density  from  the  previous 
step  into  the  current  value  of  the  solid  phase  bulk  density  (predicted  by 
Equation  (4)).  Since  the  time  steps  are  small,  and  the  changes  in  particle 
density  are  very  gradual,  the  errors  resulting  from  this  approximation  are  not 
si gni ficant. 


4.  Gas  Production  Rate 

The  burning  particles  produce  gas  according  to: 
rc  =  7^  (1  -  ♦)  Pp? 


(A-3) 


where  the  particle  surface  to  volume  ratio  for  the  spheres  is 
surface  burning  rate  was  described  as: 


3 

—  and  the 
rp 


r  = 


(A-4) 


with  b  and  n  being  known  imposed  constants 


5.  Gas  Viscosity  Coefficient 


The  viscosity  of  the  gas  phase  was  assumed  to  be  a  standard  function  of 
temperature  only,  and  was  given  as: 


u 


9 


0.65 


(A-5) 


A  relation  which  considers  the  coefficient  as  a  function  of  both  temperature 
and  pressure  would  have  been  more  desirable  here  because  the  extreme  pressures 
will  have  some  influence  on  the  viscosity. 


6,  Interphase  Drag 

The  equations  used  to  predict  the  drag  forces  between  the  gas  and  par¬ 
ticles  were: 


P 

where  the  friction  factor  f  is  defined  as: 


f 

P9 


[150  +  1.75  [y—)0'87] 


(A-6) 

(A-7) 


{ A-8 ) 


Here  Re  is  the  flow  Reynolds  number  based  upon  the  particle  diameter: 


where  V  represents  the  relative  velocity  of  the  gas  as  it  moves  over  the 
solids,  such  that: 


v  - 1  <v  up)2+  (v  wP)2]1/2 


(A-10) 


This  drag  relationship  was  derived  from  semi -empirical  results  by  Jones  and 
Krier  (Reference  24)  for  steady  flow  over  a  packed  beJ  of  soherical  heads.  It 
assumes  a  constanc  porosity  and  is  valid  for  Reynolds  number  ranging  from  733 
to  126670.  Although  the  problem  presented  in  this  study  is  unsteady  with 
variable  porosity,  the  Jones  and  Krier  correlation  produces  reasonable  values 


for  the  interphase  drag  in  the  two  dimensional  model. 


7.  Convective  Heat  Transfer 


The  rate  of  heat  transfer  between  the  gas  and  the  particles  is  described 


0  -  —  (1-*)<T  •  T  )  h 
rp  9  P  P9 


(A— 11) 


where  the  heat  transfer  coefficient  hpg  is  the  classical  Denton  formula  also 
used  by  Krier  and  Kezerle  (Reference  7): 


h  -  0.53  -S.  Re0i7  Pr0,33 
P9  rp 


(A- 12) 
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APPENDIX  B 


FINITE  DIFFERENCE  FORM  OF  THE  CONSERVATION  EQUATIONS 

Numerical  computation  methods  must  be  employed  to  solve  the  unsteady, 
two-phase  conservation  equations  derived  in  Section  III,  After  the  domain  has 
been  divided  into  discrete  nodes,  the  governing  partial  differential  equations 
are  rewritten  in  finite  difference  form  using  the  explicit,  centered  time, 
centered  space,  leapfrog  scheme.  The  resulting  algebraic  relations  are 
applied  at  all  node  points,  producing  a  system  which  can  be  solved  by  a  digi¬ 
tal  computer. 

The  following  shorthand  notation  is  used  in  the  finite  difference  equa¬ 
tions  whicn  are  presented  here.  If  “A“  is  considered  to  be  a  generalized 
dependent  variable  while  "xu  is  an  independent  variable,  then: 


A  ax  -  A  ax 

.  *  _  2  2 
6  A  55 

X  AX 


A  +  A  -  2A 

,  .  X+AX  X-AX  X 

6  it*  = - - 

*X  (AX)Z 

\+  1 5.  +  1 1 

t  x  _  c  c. 


(B-l) 

(8-2) 

(B— 3) 


Thus  the  conservation  equations  in  finite  difference  form  are: 


1.  Conservation  of  Mass 

a.  Gas  Phase  Continuity  (applied  at  scalar  rodes) 
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Particle  Phase  Continuity  (applied  at  scalar  nodes) 


5tp2t  =  "  r  rup)  "  5z^p2  wp^  ‘  rc 


(B- 


Conservation  of  Momentum 
Gas  Phase  Momentum 

(1)  Radial  Direction  (applied  at  radial  velocity  nodes) 


6Jpir  ug}  =  “  7  Mplr  uqP  “g^  “  *2^1^  “g*  wgP) 


+  pgr  <M?  Mrug^  +  Vg  +  T  Mr  Mrug>  +  5zWg» 


+  Tc  up  -  Dr  -  MV 


(B- 


(2)  Axial  Direction  (applied  at  axial  velocity  nodes) 


M’l2  "gl  ■  -  7  sr^lf  r  V  "grl  ■  'z^  *9Z  "g!> 


♦  g*  |i  «r(r«rV  *  «22».  4  «z  ir  «r(rV  *  VgU 


b.  Particle  Phase 

(1)  Radial  Direction  (applied  at  radial  velocity  nodes) 


Mp2  M  "  r  Mp2r  up  “p  )  '  sz^2  upZ  wpr) 


P  P 


P  P 


V  up  +  °r  ‘  MV1**)] 


(2)  Axial  Direction  (applied  at  axial  velocity  nodes) 


5J°2Z  wp)  "  7  srf“2^  r  upZ  "p)  ■  5zf°2wpZ  “pZ) 
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}.  Conservation  of  Energy 

i.  Gas  Phase  (applied  at  scalar  nodes) 


6t(prg)  3  '  r  6r^irugr^gr+  “  MplS(V+ 


+  Pg  l2  (^rug)2  +  +  (6zwg)2] 

-  !  !sr“g  *  V  *  V/  ‘  (‘A  *  Vj)! 
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+  Mr  Mr  ugr]  +  6z  "o'*]  +  ~p~  Mr  5r  *qZ) 
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(B-10) 


b.  Particle  Phase  (applied  at  scalar  nodes) 

F_r 


M»2y  -  7  5,  F/V<V*  ^7)1  -  «ztV«B(V  *  ^-)i 


p"  p 
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(8-11 ) 


In  Equations  (8-10)  and  (B-ll)  the  total  energies  are  defined  as: 


CVT  *  lip2.  ISp 

(8-12) 

=  c  T  r  +  (i?  +  lill2 

V  ?  1 

(8-13) 

.cT^.iii2.-2 

V  2  2 

(B-14) 

where  a  subscript  is  added  to  denote  gas  or  particle  phase.  In  addition,  the 
averages  of  the  interphase  drag  terms  must  be  described  as: 
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Ha. 

_JjL 
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<V-  V>  fP9 

(B-15) 

f»V'  V  )  fpg 

(B-16) 

while  the  remainder  of  the  source/sink  terms  present  no  significant  problems 
resuiting  from  discretization  of  the  system. 
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APPENDIX  C 


COMPUTER  CODE  LISTING 


PROGRAM  DDT2DA ( INPUT , OUTPUT , TAPE6 , TAPE7 , TAPES ) 

C 

C 

c 

C  WRITTEN  BY  LAURENCE  S.  SAMUELSON  AND  HERMAN  KRIER 

C  UNIVERSITY  OF  ILLINOIS  AT  URB ANA -CHAMPAIGN 

C  LATEST  VERSION— AUGUST  1984 
C 

c 

c 

C  THIS  IS  DDT2DA,  A  PROGRAM  WHICH  MODELS  A  DAMAGED  GENERAL  PURPOSE 

C  WARHEAD  IN  TWO  DIMENSIONS  AND  INCLUDES  THE  EFFECTS  OF  GAS 

C  VISCOSITY.  THE  PROGRAM  USES  A  CENTERED  TIME,  CENTERED  SPACE, 

C  EXPLICIT  DIFFERENCING  SCHEME  ON  A  STAGGERED  GRID.  TO  REDUCE  THE 

C  TIME  SPLITTING  OSCILLATIONS  INHERENT  IN  CENTERED  TIME  METHODS, 

C  A  SMOOTHING  STEP  IS  ALSO  INCLUDED. 

C 

C  THE  MAIN  PROGRAM  SETS  THE  SIZE  OF  THE  TIME  STEPS  AND  CONTROLS 
C  THE  FLOW  OF  INFORMATION  THROUGH  THE  SUBROUTINES. 

C 

C  SUBROUTINE  READA  READS  IN  THE  INPUT  PARAMETERS  FROM  AN 
C  EXTERIOR  DATAFILE  CALLED  DATA2D.  THESE  VALUES  ARE  THEN 
C  PRINTED  TO  AN  OUTPUT  FILE  FOR  CONFIRMATION  AND  INDEXING 
C  PURPOSES. 

C 

C  SUBROUTINE  INITIAL  LOADS  THE  VARIABLE  ARRAYS  WITH  TIE  INITIAL 
C  CONDITIONS  FOR  THE  BORN. 

C 

C  SUBROUTINE  VARSET  RESETS  AMD  UPDATES  THE  SECONDARY  VARIABLES 

C  BASED  UPON  THE  RESULTS  FROM  THE  PREVIOUS  TIME  STEP.  IT  ALSO 

C  PREPARES  THE  PRIMARY  VARIABLES  FOR  THE  NEXT  TIME  STEP. 

C 

C  SUBROUTINE  SWEEP  PERFORMS  THE  ACTUAL  FINITE  DIFFERENCE 

C  INTEGRATIONS  ON  THE  PRIMARY  VARIABLES.  ALL  THE  CALCULATIONS 

C  ARE  PERFORMED  ON  A  SINGLE  NODE  BEFORE  ADVANCING  TO  THE  NEXT 

C  POSITION.  THE  SUBROUTINE  SWEEPS  DOWN  THE  LENGTH  OF  THE  BCD 

C  BEFORE  MOVING  MOVING  TO  A  NEW  RADIAL  POSITION. 

C 

C  FUNCTION  AVG  AVERAGES  THE  VALUES  AT  THE  INPUT  NODE  AND  THE 
C  PREVIOUS  NODE  IN  THE  INPUT  DIRECTION. 

C 

C  FUNCTION  DBLAVC  AVERAGES  THE  VALUES  OF  FOUR  MODES 
C 

C  SUBROUTINE  DUMP  CONTROLS  THE  FORMAT  OF  THE  OUTPUT. 

C 

C  SUBROUTINE  DUNF2  PRINTS  RATIOS  OF  DRAG  AND  VISCOUS  STRESS  TERMS. 
C 
C 
C 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


THE  VARIABLES  USED  IN  THE  PROGRAM  ARE: 

PRIMARY  VARIABLES 

RHOI  BULK  GAS  DENSITY 

RH02  BULK  PARTICLE  DENSITY 

UG  GAS  RADIAL  VELOCITY  COMPONENT 

UP  PARTICLE  RADIAL  VELOCITY  COMPONENT 

VG  GAS  AXIAL  VELOCITY  COMPONENT 

VP  PARTICLE  AXIAL  VELOCITY  COPONENT 

EG  GAS  TOTAL  ENERGY 

EP  PARTICLE  TOTAL  ENERGY 

SECONDARY  VARIABLES 

RHOG  GAS  PHASE  DENSITY 

RHOP  PARTICLE  PHASE  DENSITY 

PG  GAS  PHASE  PRESSURE 

PP  PARTICLE  PHASE  PRESSURE 

PGPHI  BULK  GAS  PRESSURE 

PP1MF  BULK  PARTICLE  PRESSURE 

TG  GAS  TEMPERATURE 

TP  FARTICLE  PRESSURE 

PHI  POROSITY 

RADP  INSTNTAXEOUS  PARTICLE  RADIUS 

IGN  CONDITION  OF  PARTICLE  ( INTACT, BURNING, OR  BURNED  OUT) 

MUG  GAS  VISCOSITY 

GAMMAC  GAS  GENERATION  RATE  DURING  COMBUSTION 

DRAGR  RADIAL  COMPONENT  GAS-PARTICLE  INTERACTION  DRAG 

DRACZ  AXIAL  COMPONENT  GAS-PARTICLE  INTERACTION  DRAG 

QDOT  HEAT  TRANSFER  RATE  FROM  GAS  TO  PARTICLE  PHASE 

ARVISR  ARTIFICIAL  VISCOSITY  IN  THE  RADIAL  DIRECTION 

ARVISZ  ARTIFICIAL  VISCOSITY  IN  THE  AXIAL  DIRECTION 

OTHER  VARIABLES 

RLEN  ACTUAL  RADIAL  LENGTH  OF  BED 

ZLEN  ACTUAL  AXIAL  LENGTH  OF  BED 

DELR  RADIAL  DIRECTION  DIFFERENCING  SPACE  STEP 

DELZ  AXIAL  DIRECTION  DIFFERENCING  SPACE  STEP 

DELT  DIFFERENCING  TIME  STEP 

STEP  INCREMENTAL  TIME  STEP  FACTOR 

EPS  TIME  SMOOTHER  COEFFICIENT 

ARV  ARTIFICIAL  VISCOSITY  COEFFICIENT 

TOL  TRUNCATION  ERROR  TOLERANCE 

R  ACTUAL  RADIAL  POSITION 

TIME  ACTUAL  TINE 

CVG  GAS  SPECIFIC  HEAT  AT  CONSTANT  VOLUME 

CVP  PARTICLE  SPECIFIC  HEAT  AT  CONSTANT  VOLUME 

RO  GAS  CONSTANT 

MV  MOLECULAR  HEIGHT  OF  GAS 

B1  NON-IDEAL  EQUATION  OF  STATE  CONSTANT 

TGO  INITIAL  BULK  GAS  TEMPERATURE 

TPC  INITIAL  BULK  PARTICLE  TEMPERATURE 

TCH  INITIAL  LOCAL  MAXIMUM  TEMPERATURE 

TIGN  TEMPERATURE  NECESSARY  TO  IGNITE  PARTICLES 

B  BURNING  RATE  PROPORTIONALITY  CONSTANT 

BUT  BURNING  RATE  INDEX 

EPT  PROPELLANT  ENERGY 

RHOPO  INITIAL  PARTICLE  DEK3ITY 
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PHIO  INTTIAL  POROSITY 

PGO  INITIAL  BULK  GAS  PRESSURE 

PCH  INITIAL  LOCAL  MAXIMUM  PRESSURE 

PWALL  MAXIMUM  PRESSURE  AT  THE  HALLS 

RADPO  INITIAL  PARTICLE  RADIUS 

KO  PARTICLE  BULK  MODULUS 

HUGO  INITIAL  GAS  VISCOSITY 

PR  PRANDTL  NUMBER 

INIT  DETERMINES  THE  TYPE  OF  INITIATION 

(ONE  OR  TWO  DIMENSIONAL) 

TOGGLE  TOGGLES  GAS  VISCOSITY  EFFECTS  ON  AND  OFF 
ENDR  BOUNDARY  NODE  IN  THE  RADIAL  DIRECTION 

ENDZ  BOUNDARY  NODE  IN  THE  AXIAL  DIRECTION 

STOP IT  WARNING  FLAG  FOR  IMPOSSIBLE  CONDITIONS  OCCUR  I NG 
COUNT  INDICATES  THE  NUMBER  OF  SWEEPS  THROUGH  THE  BED 
MAX  MAXIMUM  NUMBER  OF  INTEGRATION  SWEEPS  PERMITTED 
RPRINT  DETERMINES  THE  NUMBER  OF  RADIAL  NODES  PRINTED 
ZPRINT  DETERMINES  THE  NUMBER  OF  AXIAL  NODES  PRINTED 
TPRINT  DETERMINES  THE  NUMBER  OF  TIME  STEPS  BETWEEN  PRINTING 
NTGR-  NUMBER  OF  NODES  AT  INITIAL  TEMPERATURES/PRESSURES 
NPGZ  ABOVE  THE  INITIAL  BULK  TEMPERATURE/PRESSURE 
LOCAL  VARIABLES 

A  SPEED  OF  SOUND 

BURNED  AVERAGE  CONDITION  OF  PARTICLES  AT  END  OF  BED 
PHIM  1-POROSITY 

ROOT  BURNING  RATE 

V  TOTAL  GAS-PARTICLE  RELATIVE  VELOCITY 

RE  REYNOLDS  NUMBER 

FPG  DRAG  COEFFICIENT 

KG  GAS  THERMAL  CONDUCTIVITY 

HTC  GAS  TO  PARTICLE  HEAT  TRANSFER  COEFFICIENT 

CGI-  INTERMEDIATE  QUANTITIES  IN 

EP6  CONTINUITY,  MOMENTUM,  AND  ENERGY  EQUATIONS 

UGRH01-  INTERMEDIATE  QUANTITIES  IN 
EPRH02  TIME  SMOOTHER  EQUATIONS 
ECCHEH  GAS  CHEMICAL  ENERGY  RELEASED  DURING  COMBUSTION 
EPCHEM  PARTICLE  CHEMICAL  ENERGY  RELEASED  DURING  COMBUSTION 
RLOC  ACTUAL  R  LOCATION  FOR  PRINTING 

ZLOC  ACTUAL  Z  LOCATION  FOR  PRINTING 

PGSI  GAS  PRESSURE  IN  GPA 

RADPSI  PARTICLE  RADIUS  IN  MICRO- METERS 

TIMOUT  TIME  OF  PRINTING  IN  M  iCRO-SECONDS 
DTOUT  SIZE  OF  TIME  STEP  IN  MICRO-SECONDS 
PVALOT  MAXIMUM  WALL  PRESSURE  IN  GPA 
I.J.K  DO  LOOP  COUNTERS 


**  LOCAL  VARIABLES  ** 
REAL  A 

INTEGER  I,  J,  BURNED 


**  COMMON  VARIABLES  M 
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c 

REAL  RH01(0:6,0:61 ,3) »  RH02(0:6,0:6l ,3) ,  UG(0:6, 0:61.3), 

&  UP(0:6, 0:61,3),  HG<0 : 6 ,0 :6l , 3) #  VP(Q:6,0:6l ,3), 

&  EG(0:6, 0:61,3),  EP(0:6 ,0:61 , 3) 

REAL  RH0G(0:6,0:61),  RH0P{0:6,0:61),  PG(0:6,0:61),  PP(0:6,0:6l), 
A  PGPHI(0:6,0;6l),  PP1«P(0:6,0:6l ) ,  TG(0:6,0:6l ) ,  TP(0:6,0:6l ), 

&  PHI(0:6 ,0:61 ) ,  RADP(0:6,0:6l ) ,  HUG(0:6,0:6l), 

&  GAMMAC(0:6,0:6l) ,  DRAGR(0:6,0:6l ) ,  DRAG2(0:6,0:6l ), 

A  QDOT(0:6,0:61),  ARVISR(0:6,0:6l),  ARVJSZ(0:6,Q:6l) 

INTEGER  ICN(0:6,0:61) 

REAL  RLEN,  ZLEN,  DELR,  DELZ,  DELT,  STEP,  EPS,  ARV,  TOL,  fl,  TIME, 
A  CVG,  CVP,  RO,  MW,  B1,  TGO,  TPO,  TCH,  TIGN,  B,  ENT,  EPT,  RHOPO, 

A  PHIO,  PGO,  PCH,  PWALL,  RADPO,  KO,  HUGO,  PR 
INTEGER  INIT,  TOGGLE,  ENDR,  ENDZ,  STOPIT,  COUNT,  MAI, 

A  RPRINT,  ZPRINT,  TPRINT,  NTGR,  NTGZ,  NTPR,  NTPZ,  NPGR,  NPGZ 
COMMON  /PRI/  RH01 ,  RH02,  UG,  UP,  WG,  WP,  EG,  EP 
COMION  /SEC/  RHOG,  RHOP,  PG,  PP,  PCPHI,  PP1MF ,  TG,  TP,  PHI, 

A  RADP,  IGH,  MUG,  GAMMAC,  DRAGR,  DRAGZ,  QDOT,  ARVISR,  ARVISZ 
COMMON  //  RLEN,  ZLEN,  DELR,  DELZ,  DELT,  STEP,  EPS,  ARV,  TOL,  R, 

A  TIME,  CVG,  CVP,  RO,  MW,  B1,  TGO,  TPO,  TCH,  TIGN,  B,  BNT,  EPT, 

A  RHOPO,  PHIO,  PGO,  PCH,  PWALL,  RADPO,  KO,  HUGO,  PR,  INIT,  TOGGLE, 
A  ENDR,  ENDZ,  STOPIT,  COUNT,  MAI, 

A  RPR I NT,  ZPRINT,  TPRINT,  NTGR,  NTGZ,  NTPR,  NTPZ,  NPGR,  NPGZ 

C 

C 

CALL  READA 
CALL  INITIAL 
TIME=0.0 
COUNT=Q 
DELTsl.OE-8 
CALL  DUMP 
CALL  VARSET 
DO  10  COUNT: 1, MAI 
CALL  SWEEP 
C 

C  M  CHECK  FOR  BURNED  OUT  CONDITIONS  ** 

C 

BURMEDsO 
DO  20  1:1, ENDR 
BURNED: 3URKED+ IGM ( I , ENDZ- 3 > 

20  CONTINUE 

IF  ( <BURMED/ENDR)ffENDR.G£. 1 )  THEN 
CALL  DUMP 
STOP 
END  IF 
C 

IF  ( (CQUNT/TPRIMT)*TPRINT.£Q.CGUNT }  CALL  DUMP 
CALL  VARSET 
?IHE=?IME+DELT 
C 

C  **  COMPUTE  MEW  TIME  STEP  ** 

C 

DELT=5.0E-7 
DO  30  1=1, ENDR 
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-*v* 

*-V; 


DO  40  J=1,ENDZ 

A=SQRT(PG(I,J)/RHOG( I,J)*( ( 1.Q+2.Q*B1*RH0G(I,J)) 

A  /( 1 .0+B1*RH0G( I , J) )+R0/MW/CVG*( 1 .0+B1«RHOC< I ,  J) ) ) ) 

A=SQRT( K0/RH02 ( I , J , 2 ) ) 

DTR=STEP*D£LR/A 
DTZ=STEP*DELZ/A 
IF  (OTR.LT.DELT)  DELT-DTR 
IF  (DTZ.Lt.DELT)  DELTsDTZ 
CONTINUE 
CONTINUE 
CONTINUE 

»»  IF  THIS  POINT  IS  REACHED,  THE  PROGRAM  HAS  EXCEEDED  ** 

«  THE  MAXIMUM  NUMBER  OF  TIME  STEPS  ALLOWED  M 

PRINT  100 

FORMAT!/'  THE  MAXIMUM  NUMBER  OF  TIME  STEPS  HAS  BEEN  REACHED' ) 

STOP 

END 


SUBROUTINE  REA DA 

**  COMMON  VARIABLES  “ 

REAL  RHO1(O:6,0j61,3).  RH02(0:S,0:6l,3),  UG{0:6, 0:61,3), 

&  UP(0: 6,0:61 ,3) ,  WG(0:6,0:6l ,3) ,  WP(0:6,0:6l , 3) , 

&  EG(0: 6,0:61, 3),  EP(0:6,0:6l ,3) 

REAL  RHOG(0:6,O:6l),  RH0P(0:6,0:6l),  PG(0:6,0:6l),  PP(0:6,0:6l), 
&  PGPHI (0:6,0:61 ) ,  PPIlff (0:6,0:61 ) ,  TG(0:6,0:6l),  TP(0:6,0:Sl), 

&  PHI(0:6,0:61 ) ,  RADP(0:S,0:61) ,  MUG{0:6,0:61), 

A  GAMMAC(0:6 ,0:61 ) ,  DRAGR(0:6,0:61 ),  DRAGZ(0:6,0:6l ) , 

A  QD07(C:6 ,0:61 ) ,  ARVISR(0:6,0:6l),  AJWISZ(0:6,0:6l } 

INTEGER  ICN(0:6,0:6l) 

REAL  RLEN,  ZLEN,  DELK,  DELZ,  DELT,  STEP,  EPS,  ARV,  TOL,  R,  TIME, 
A  CVG,  CVP,  RQ,  MW,  B1,  TOO,  TPC,  TCH,  TIOI,  B,  BUT,  EPT,  RHOPO, 

A  PH 10,  POO,  PCH,  PVALL,  RADPO,  10,  HUGO,  PR 
INTEGER  INIT,  TOGGLE,  ENDR,  ENDZ,  STOPIT,  COUNT,  MAX, 

A  RPR I NT,  ZPRINT,  TPRINT,  NTGR,  NTGZ,  NTPR,  NTPZ,  NPGR,  NPCZ 
COWCN  /PRI/  RH01,  RHD2,  UG,  UP,  NG,  HP,  EG,  EP 
COMMON  /SEC/  RHOC,  RHOP,  PG,  PP,  PGPHI,  PP1HF,  TG,  TP,  PHI, 

A  RADP,  IGN,  MUG,  GAMMAC ,  DRACH,  DRAGZ,  QDOT,  ARV I SR,  ARVISZ 
COMHQN  //  RLEN,  ZLEN,  DELR,  DELZ,  DELT,  STEP,  EPS,  ARV,  TOL,  R, 

A  TIME,  CVG,  CVP,  HO,  MU,  B1 ,  TGO,  TPO,  TCH,  TIOI,  B,  BUT,  EPT, 

A  RHOPO,  PHIO,  PGO,  PCH,  PHALL,  RADPO,  KO,  HUGO,  PR,  INIT,  TOGGLE, 
A  ENDR,  ENDZ,  STOPIT,  COUNT,  MAX, 

A  RPRINT,  ZPRIHT,  TPRINT,  NTGR,  NTGZ,  NTPR,  NTPZ,  UPGR,  MPGZ 


READ  1,  RLEN 
READ  2,  ZLEN 
READ  3,  ENDR 
READ  4,  ENDZ 
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READ  5,  STEP 
READ  6,  EPS 
READ  7,  ARV 
READ  8,  TOL 
READ  9,  MAX 
READ  10,  CVG 
READ  11,  CVP 
READ  12,  RO 
READ  13,  MV 
READ  14,  B1 
READ  15,  TGO 
READ  16,  TPO 
READ  17,  TCH 
READ  18,  NTGR 
READ  19,  MTGZ 
READ  20,  MTPR 
READ  21,  NTPZ 
READ  22,  PGO 
READ  23,  PCH 
READ  24,  VPGR 
READ  25,  MPGZ 
READ  26,  TIC3I 
READ  27,  S 
READ  28,  BNT 
READ  29,  KO 
READ  30,  EPT 
READ  31.  RHOPO 
READ  32,  PHIO 
READ  33,  RADPQ 
READ  34,  HUGO 
READ  35,  PR 
READ  36,  TP R I 47 
READ  37,  RFRIWT 
READ  38,  IPRI4T 
READ  39,  IMIT 
READ  40,  TOGGLE 

1  FORMAT* //19X,G33.l6) 

2  FORMAT* 18X.G33-16) 

3  FORHAT(/ 191,13) 

4  FORMAT* 18X. 13) 

5  FOfWAT(29X,G33-l6) 

6  FORMAT (2YX,G33.l6) 

7  FORMAT* 33X.C33.l6) 

8  F0RMAT(36X,C33-16) 

9  FORMAT* 37X, 15) 

10  FORMAT* /28X.G33. 16) 

11  FORMAT* 33X.G33- 16) 

12  FORMAT* 35X,G33-l6) 

13  FORMAT* 29X.G33* 16) 

14  FORMAT* 44X.C33- 16> 

15  FORMAT (/33X.G33- ^6) 

16  FORMAT* 38X, 033.16) 

17  F0HSAT*38X,G33.16) 

18  FORMAT* /19X, 13) 
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19  FORMAT* 19* *13) 

20  FORMAT* /19X, 13) 

21  FORMAT* 18X, 13) 

22  F0RMAT*34X,G33.16) 

23  FORMAT(43X,G33.l6> 

24  FORMAT* /19X, 13) 

25  FORMAT* l8Xf 13) 

26  FO»«AT(/46X,G33.t6) 

27  FORMAT*38XtG33-l6) 

28  FORMAT* 19X, 033.16) 

29  FORMAT* 34X,G33- 16) 

30  FORMAT* 24X,G33- 16) 

31  FORMAT *32X,G33. 16) 

32  FORMAT* 17X.G33. 16) 

33  FORMAT (29X.G33. 16) 

34  FORMAT* 30X.G33. 16) 

35  FORMAT* 19X.G33. 16) 

36  FORMAT* /36X, 13) 

37  FORMAT*  38X,  13) 

38  FORMAT* 37X, 13) 

39  F0RMAT(//41X,I1) 

40  FORMAT* 45X, II) 


PRINT  100 
PRINT  101,  RLEN 
PRINT  102,  ZLEN 
PRINT  103,  EMDR 
PRINT  104,  ENDZ 
PRINT  105,  STEP 
PRINT  106,  EPS 
PRINT  107,  ARV 
PRINT  108,  TOL 
PRINT  109,  MAX 
PRINT  110,  CVG 
PRINT  111,  CVP 
PRINT  112,  RO 
PRINT  113,  MW 
PRINT  114,  B1 
PRINT  115,  TOO 
PRINT  116,  TPO 
PRINT  117,  TCH 
PRINT  118,  NTGR 
PRINT  119,  NTGZ 
PRINT  120,  NTPR 
PRINT  121,  NTPZ 
PRINT  122,  POO 
PRINT  123,  PCH 
PRINT  124,  WPGR 
PRINT  125,  NPGZ 
PRINT  126,  TIGN 
PRINT  127,  B 
PRINT  128,  BNT 
PRINT  129,  XO 
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PRINT  130,  EPT 
PRINT  131,  RHOPO 
PRINT  132,  PHIO 
PRINT  133,  RADPO 
PRINT  134,  HUGO 
PRINT  135,  PR 
PRINT  136,  TPRIMT 
PRINT  137,  RPRINT 
PRINT  1 38 ,  SPRINT 
IF{ INIT.EQ. 1 )  PRINT  139 
IF(INIT.EQ.O)  PRINT  239 
IF(TOGCLE.EQ.I)  PRINT  140 
IF ( TOGGLE. EQ.O)  PRINT  240 

100  FORMAT!  'IV  INPUT  PARAMETERS  FOR  DDT2DA ' ) 

101  FORMAT! 'O', 'RADIAL  LENGTH  !CH)=' ,G12.4) 

102  FORMAT!'  ', 'AXIAL  LENGTH  <CM )  =  ’ .G12.4) 

103  FORMA?!'  ', 'NUMBER  OF  SCALAR  NODES--'/ 

A  '  RADIAL  DIRECTIONS  M3) 

104  FORMAT!'  AXIAL  DIRECTION* ' , 13) 

105  FORMAT!'  '.'INCREMENTAL  TIffi  STEP  FACTORS', G12. 4) 

106  FORMAT! '  ','TIKE  SMOOTHING  COEFFICIENT* ' ,G12. 4) 

107  FORMAT!'  ', 'ARTIFICIAL  VISCOSITY  COEFFICIENT*' ,G12. 4) 

108  FORMAT!'  'TRUNCATION  ERROR  TOLERANCE  (CM/SEC)* • ,G12. 4) 

109  FORMAT!'  MAXIMUM  NUMBER  OF  TIME  STEPS  ALLOWED* *,I5) 

110  FORMAT! 'O', 'GAS  SPECIFIC  HEAT  (ERG/G/K)*' ,G12.4) 

111  FORMAT!'  ', 'PARTICLE  SPECIFIC  HEAT  (ERG/G/K)* ' ,G12. 4) 

112  FORMAT!'  ', 'UNIVERSAL  GAS  CONSTANT  (ERG/MOL/K)* ' ,G12. 4) 

113  FORMAT!'  ', 'GAS  MOLECULAR  NEIGlfr  {G/M0L)=',G12. 4) 

114  FORMAT!'  ', 'NON-IDEAL  EQUATION  OF  STATE  CONSTANT  (CC/G)s * ,G12.4> 

115  FORMAT! 'O’, 'INITIAL  BULK  GAS  TEMPERATURE  (K)s',G12.4) 

116  FORMAT!'  '.'INITIAL  BULK  PARTICLE  TEMPERATURE  (K)=*,G12.4) 

117  FORMAT!'  MAXIMUM  LOCAL  INITIAL  TEMPERATURE  (I)s\G12.4) 

118  FORMAT!'  ' .  'NUMBER  OF  SCALAR  NODES  ABOVE  BULK  GAS  TEMPERATURE—'/ 
A  '  V  RADIAL  DIRECTION* ’,13) 

119  FORMAT! '  ','  AXIAL  DIRECTION* M3) 

120  FORMAT! '  ’ , 

A  'NUMBER  OF  SCALAR  NODES  ABOVE  BULK  PARTICLE  TEMPERATURE—'/ 

4  '  ',’  RADIAL  DIRECTION* *,I3) 

121  FORMAT!’  AXIAL  DIRECTION*' , 13) 

122  FORMAT!'  ',' INITIAL  BULK  PRESSURE  (DYN/CM/OO* \G12.«> 

123  FORMAT!'  '.'MAXIMUM  LOCAL  INITIAL  PRESSURE  !DXM/CM/CM)= • .G12.4) 

124  FORMAT!’  ', 'NUMBER  OF  SCALAR  NODES  ABOVE  BULK  PRESSURE—'/ 

A  '  RADIAL  DIRECTION* M3) 

125  FORMAT!'  V  AXIAL  DIRECTION*' .13) 

126  FORMAT! '0',' TEMPERATURE  NECESSARY  TO  IGNITE  PARTICLES  (K)=',G12.4) 

127  FORMAT!'  ’.'BURNING  RATE  PROPORTIONALITY  CONSTANT* ' ,G12.4) 

128  FORMAT!’  BURNING  RATE  INDEX* ’ ,G12. 4) 

129  FORMAT!'  ' , ’PARTICLE  BULK  MODULUS  (DYN/QC/CM)* ' ,G12.4) 

130  FORMAT!'  ' , 'PROPELLANT  ENERGY  (ERG)* ' ,G 12. 4) 

131  FORMAT! '  '.'INITIAL  PARTICLE  DENSITY  (G/CC)* • ,G12.4) 

132  FORMAT! '  'INITIAL  POROSITY* ' ,G12. 4) 

133  FORMAT!'  '.'INITIAL  PARTICLE  RADIUS  (CM)*’ ,G12.4) 

134  FORMAT!'  'INITIAL  GAS  VISCOSITY  {POISE}*' ,G12.4) 

135  FORMAT!'  ','GAS  PRANDTL  NUMBER* ' ,C12. 4) 


136  FORMAT COVHUMBER  OF  TIME  STEPS  BETWEEN  OUTPUT=\I3) 

137  FORMAT  ( '  NUMBER  OF  RADIAL  MODES  BETWEEN  OUTPUT:  M3) 

138  FORMAT  ( '  V NUMBER  OF  AXIAL  NODES  BETWEEN  OUTPUT:  M3) 

139  FORMATC -*,'»*  INITIATION  TAKES  PLACE  TWO-DIMENSIONALLY  **') 

239  FORMAT{ •-',«**  INITIATION  TAKES  PLACE  ONE-DIMENSIONALLY  **' ) 

140  FORMATC VISCOUS  CAS  EFFECTS  ARE  INCLUDED  **> ) 

240  FORMATC O', '•*  VISCOUS  GAS  EFFECTS  ARE  NOT  INCLUDED  **• ) 
RETURN 

END 


C 

C 

c 

c 

c 

c 

c 

c 

c 


c 

c 


c 

c 

c 


SUBROUTINE  INITIAL 
»*  LOCAL  VARIABLES  H 
INTEGER  I,  J,  K 
«  COMMON  VARIABLES  ** 


REAL  RH01(0:6,0:61,3),  RH02(0:6,0:6l ,3),  UG(0:6,0:6l,3), 

&  UP(0:6,0:6l ,3) ,  WG(0:6,0:61 ,3) ,  WP(0:6,0:6l,3), 

&  EG(0:6,0:61,3),  EP(0:6,0:61 ,3) 

REAL  RHOG(0:6,0:61 ) ,  RH0P(0:6,0:6l),  PG(0:6,0:61 ) ,  PP(0:6,0:6l ), 
&  PGPHI(0:6 ,0:61 ) ,  PP1MF(0:6,0:6l ) ,  TG(0:6,0:6l),  TP(0:6,0:6l) , 

A  PHI{0:6,0:61 ) ,  RADP(0:6,0:6l ) ,  MUG(0:6,0:61 ), 

&  GAHMAC(0:6,0:61),  DRAGR(0:6,0:61 ) ,  DRAGZ(0:6,0:6l ) , 

&  QDOT(0:6,0:61),  ARVISR{0:6,0:61 } ,  ARVISZ(0:6,0:6l ) 

INTEGER  IGN(0:6,0:61) 

REAL  RLEN,  ZLEN,  DELR,  DELZ,  DELT,  STEP,  EPS,  ARV,  TOL,  R,  TIME, 
4  CVG,  CVP,  RO,  MW,  B1,  TOO,  TPO,  TCH,  TIGN,  B,  BNT,  EPT,  RHOPO, 

A  PHIO,  PGO,  PCH,  PWALL,  RADPO,  KO,  HUGO,  PR 
INTEGER  INIT,  TOGGLE,  ENDR,  ENDZ,  STOPIT,  COUNT,  MAX, 

A  RPRINT,  ZPRINT,  TPRINT,  NTGR,  NTGZ,  NTPR,  NTPZ,  MPGR,  NPGZ 
COMWN  /PRI/  RH01 ,  RH02,  UG,  UP,  WG,  WP,  0),  EP 
COMMON  /SEC/  RHOG,  RHOP,  PG,  PP,  PGPHI,  PP1MF,  TG,  TP,  PHI, 

A  RADP,  IGN,  MUG,  GAMfAC,  DRAGR,  DRAGZ,  QDOT,  ARVISR,  ARVISZ 
COMMON  //  RLEN,  ZLEN,  DELR,  DELZ,  DELT,  STEP,  EPS,  ARV,  TOL,  R, 

4  TIME,  CVG,  CVP,  RO,  MW,  B1,  TGO,  TPO,  TCH,  TIGN,  B,  BNT,  EPT, 

A  RHOPO,  PHIO,  PGO,  PCH,  PWALL,  RADPO,  KO,  HUGO,  PR,  INIT,  TOGGLE, 
A  ENDR,  ENDZ,  STOPIT,  COUNT,  MAX, 

A  RPRINT,  ZPRINT,  TPRINT,  NTGR,  NTGZ,  NTPR,  NTPZ,  NPCR,  NPGZ 


EXTERNAL  AVG 
STOPIT=0 
QELR=RLEN/ENDR 
DELZ=ZLEN/ENDZ 
DO  10  I=0,ENDR+1 
DO  20  J=0,ENDZ+1 

«  POROSITY  H 


PHI(I,J)=PHIO 


87 


”  TEMPERATURE  -- 
TG(I,J)=TGO 

IF(IKIT.EQ. 1 .AWD.(I.LE.MTGR.AND. J.LE.HTGZ) )  THEN 
IF(J+(1.0*NTGZ)/NTGR*I.LE  NTGZ) 

4  TG ( I , J ) = ( NTGZ*NTGR*TCH-NTGZ# ( TCH-TGO ) * ( I- 1 ) -NTGR 

4  *(TCH-TGO)*( J-1 ) )/(NTGZ*NTGR) 

TG( I ,J)=TGO+ (TCH-TGO )*( (NTGR*1 .0-( 1-1 ) )/NTGR*1 .0)##14.0 
4  *((NTGZ*1.0-(J-1))/NTGZ*1.0)*»iH.o 

ELSE 

IF( J.LE.NTGZ) 

4  TG{ I ,J)=7CH- (TCH-TGO )*(J- 1 )/NTGZ 

4  TGd,J)=TGO+(TCH-TG0)*((NTGZ*1.0-(J-1))/NTGZ*1.0)**l4.0 

END  IF 
TP( I , J)=TPO 

IF(INIT.EQ. 1 . AND.( I .LE.NTPR.AND. J.LE.NTPZ) )  THEN 
IF(J+(1.0*NTPZ)/NTPR*I.LE.NTPZ) 

4  TP(I,J)=(NTPZ*NTPR*TCH-NTPZ*(TCH-TPO)*(I-1)-NTPR 

4  #(TCH-TPO)*( J-1 ) )/(NTPZ*NTPR) 

TP( I,J)=TPO+(TCH-TPO)*( (NTPR*1 .0-(I-1 ))/NTPR*1 .0)**14.0 
4  *( (NTPZ*1 .0-( J-1 ) )/NTPZ*1 ,0)**14.0 

ELSE 

IF( J.LE.NTPZ) 

4  TP( I , J)=TCH-(TCH-TPO)*( J-1 )/NTPZ 

4  TP(I, J)=TPO+(TCH-TPO)*( (NTPZ*1 .0-( J-1 ))/MTPZ*1 .0)**14.0 

END  IF 

«  PRESSURE  ** 

PG(I,J)=PGO 

IF(INIT.EQ. 1 .AND. ( I.LE.NPGR.AND. J.LE.NPGZ) )  THEM 
IF( J+( 1 .0*NPGZ)/NPGR»I .LE.NPGZ) 

4  PG(  I ,  J)  =  (NPGZ*NPGR*PCH-NPGZ*(PCH-PGO)*d-1  )-NPGR 

4  »(PCH-PGO)*(J-l))/(NPGZ*NPGR) 

PG(I,J)=PG0+(PCH-PGO)*((NPGR*1.0-(I-1))/NPCR*1.0)«l4.0 
4  *( (NPGZ*1 .Q-( J-1 ) )/NPGZ*1 .0)##14.0 

ELSE 

IF( J.LE.NPGZ) 

4  PG( I , J)=PCH-(PCH-PGO)*( J-1 )/NPGZ 

4  PG(I,J)=PGO+(PCH-PGO)*((NPGZ*1.0-(J-1))/NPGZ»1.0)**14.0 

END  IF 

PP(I,J)=PG(I,J) 

PGPHId , J)=PG(  I , J)*PHI(I ,  J) 

ppiMFd,j)=pp(i,j)»(i.o-pHid,j>) 

P«ALL=PCH 
**  DENSITY  ** 

RHOP(I, J)=RHOPO 

RHOG ( I , J ) = { PC (I , J ) *HW ) / ( RQ*TG (I , J ) ) 

”  PARTICLE  SIZE  AND  GAS  PRODUCTION  «* 


A 
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>£■ 

1 


PS? 


KV 


C  >»-• 


»v 


■5$ 

1‘a1 


c 

r 


C 

C 


c 

c 

c 


c 

c 

c 


c 

c 

c 


RADP(I,J)=RADPO 
IGN(I,J)=0 
GAHHAC(I, J)=O.Q 


"  DRAG,  VISCOSITY,  AMD  HEAT  TRANSFER  RATE  « 


DRAGR( I , J)=0.0 
DRAGZ( I , J)=0.0 
MUG(I,J)=0.0 
QDOT( I , J)=0.0 


**  PRIMARY  VARIABLES  ** 


30 

20 

10 


DO  30  Ks1,3 

RHOUl,  J,K)=RROG(I,J)*PHI(I,J) 

RHO2(I,J,K)=RHOP(I,J}*(1.0-PHI(I,J)} 

UG(I,J,K)=0.0 

UP(I,J,K)=0.0 

WG( I , J,K)=0.0 

WP( I , J,K)=0.0 

EG{ I , J,K)=(TG( I, J)-TG0)*CVG 
EP(I,J,K)=(TP(I,J)-TPO)*CVP 
CONTINUE 
CONTINUE 
CONTINUE 
RETURN 
END 


SUBROUTINE  VARSET 


**  LOCAL  VARIABLES  ** 


REAL  PHIM,  ROOT,  V,  RB,  FPC-,  KG,  HTC 
INTEGER  I,  J,  K 


COMMON  VARIABLES  ** 


REAL  RH01(0:6,0:61,3),  RH02(0:6,0:61 ,3) ,  UG(0:6, 0:61,3), 

&  UP(0:6, 0:61,3),  WG(0:6,0:61 ,3),  HP(0:6,0:6l,3), 

A  EG(0:6,0:6l ,3) ,  EP(0:6,0:61 , 3) 

REAL  RHOG(0:6,0:61),  RH0P(Q:6,0:6l),  PG(C:C,0:6l) ,  PP(0:6,0:6l), 
4  PGPHI(0:6,Q:61),  PP1MF(Q:6,0:6l),  TG(0:6,G:6l),  TP(0:6,0:61), 

A  PHI (0:6 ,0:6 1 ) ,  RADP(0:6,0:6l ) ,  MUG(0:6,0:61 ) , 

&  GAKMAC(0: 6,0:61 ),  DRAGR(0:6,0:61 ) ,  DRAGZ(0:6,0:6l ), 

4  QDOT(0:6,0:61),  ARVISR(0:6,0:6l ) ,  ARVISZ(0:6,0:61 ) 

INTEGER  IGN (0:6,0:61 } 

REAL  RLEN,  ZLEM,  DELR,  DELZ,  DELT,  STEP,  EPS,  ARV,  TOL,  R,  TIME, 
&  CVG,  CVP,  RO,  MW,  B1,  TOO,  TPO,  TCH,  TIGN,  B,  BUT,  EPT,  RHOPO, 

&  PHIO,  PGO,  PCH,  PWALL,  RADPO,  KO,  HUGO,  PR 
INTEGER  INIT,  TOGGLE,  ENDR,  ENDZ,  STOPIT,  COUNT,  MAX, 

&  RPRINT,  ZPRINT,  TPRINT,  NTGR,  NTGZ,  NTPR,  NTPZ,  NPGR,  NPGZ 
COMMON  /PRI/  RH01 ,  RK02,  UG,  UP,  WG,  WP,  EG,  EP 
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0 
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c 

c 

c 


c 

c 

c 


c 

c 

c 


p 


c 

c 


0 

c 

0 


,% 
A 
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.Wi 

X;Cs 
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cotami  /sec/  rhqg,  shop,  po.  pp,  pgphx,  primp,  to,  to,  phi, 

*  RAPP,  IOI,  HUG,  OAJtOC,  DRAQR,  DRAGZ,  900T,  ARVISR,  ARIISS 
COHMOH  //  RLEE,  ZLEE,  DSLB,  PEW,  PELT,  STEP,  IPS,  ARP,  TOP,  R, 

A  Hie,  CPO,  CPP,  80,  ME,  81,  TOO,  TOO,  TOR,  TlOf,  8,  Wf,  EPT, 

*  RHPPO,  PHIO,  POO,  PCH,  PEALL,  BAPR9,  SO,  HUGO,  PR,  IEXT,  TOGGLE, 
A  BMDR,  EE3?,  STOP  IT,  COVE?,  MAS, 

A  RP8IHT,  ZPBJET ,  TPRIHT,  HTGB,  PTOZ,  RTF*,  MTOZ,  MPGR,  MfGZ 


EITEHEAL  APQ 
DO  10  I=Q,EEDR 
DO  20  J=0,EHDZ 

*•  POROSITY  APPRO# XMATIOR  *f 

PHIM=BH02 ( I , 4 , 1 ) /BHDP( I , J ) 
PHI(I,J)-1.0-PHIH 


«  TEMPERATURE  ” 

T0(I,4)=(BQ(  1 ,4, 1  )-A?G(UO,  1+1  ,4,  1 , 1  )##?,Q/2,9 
A  -AVQtHO.I^+I.I.ajtfR^/RrOI/qpStTO© 

TP(I,4)=(BP(I,4, D“AVQ(UP,X*1 f 4,1 , 1)t*gtO/#(0 
A  -APG(W,X,4+1f1,2)W2.Q/a.Q)/CPP*.TO9 

«  EOUATIORS  Op  STATE  ft 


FJIPGI 1 ,4  )=RJ»1  ( X ,  4 , 1  5/PHX  ( 1 ,4 ) 

PG{  I,J)=(RQ/M¥)tTG(If4)*RHPQ(X,4)*( 1,0+B1*RHQG(I,J)) 
PG?ai(If4)=Pq(I,4)*PHX(I,4) 

?pa,4)=PG(X,4) 

PP1MF( I l4)=PP{I,4)*PHIM 

RHOP(  1 ,4)=RlfOpO#(3?Qf(pp(  1 ,4M ,0iS)/1»*1 ,0)##{  1 .0/3.0) 

«  CORTAIHMEHT  WAX4,  PRESSURE  ft 

IF(PG(EEPE,4)tOT.PEALL>  PEAI4,=P9(1PB»,4> 

!F(PG(X,EMOZ) -G?. PEALL)  PEALL=PG(I,EEUZ> 

«  GAS  PBOPUCTXOH  ft 


IP(PHI(lt4)rGTr0,9«S)  THEE 
IQE(X,4)=? 

GANMA0Ui4)=Q*Q 

IP(TO(X,4)?0TtTG(|,4))  TO(J,4)?T0U,4) 

ELSE  XP(TOU, 4),GE,TXGE, OR, X9E(X,4).IQ,1)  THEE 
IGE(I.4)-1 

RD0T=(Bt(fYUX,4)*1  tHiE^§)t*BET)*i.|4 
RAUP(  X  ,4)=RA0P(  X  ,4>*Rm»0ffLT 
GAIWAC(li4)=3-Q/RAPPU,4)fPHIH*R80P(I,4)*ROOT 

EL&E 

IGE( 1,45=0 
QAMWtC(l,4)=9,Q 

EEP  XP 
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«  GAS  VISCOSITY  »* 

MUG ( I , J ) = HUGO* ( ( TG ( I , J ) /TGO ) **0 . 6 5 ) 

**  REYNOLDS  NUMBER  •» 

V=SQRT( (UG( I,J,1)-UP(I,J,1) )**2.0+(WG(I , J, 1 )-WP( I , J, 1 ) )**2.0) 
RE=RHO H I , J , 1 )»V«2 . 0»RADP( I , J )/MUG( I , J ) 


FPG=PHIM«*2.0/PHI(I,J)M2.0*(150.0+3.89»(RE/PHIM)«*0.87) 

DRAGR(I,J)=MUG(I,J)/(4.0*RADP(I,J)**2.0)*FPG 

*(UG(I,J,2)-UP(I,J,2)) 

DRAGZ( I , J } =MUG( I , J )/(4.0*RADP( I , J)**2.0)*FPG 
*(MG( I , J,2)-WP( I , J ,2) ) 

**  HEAT  TRANSFER  RATE  ** 

KG=MUG( I, J)*( (RO/MW)+CVG)/PR 

HTC=0 .58«KG/RADP( I , J )»(RE**0. 7 )»(PR**0. 3333) 

IF( IGN( I , J) .GT.O)  THEN 
QD0T(I,J)=0.0 
ELSE 

QDOT(I,J)=3.0«PHIM/RADP(I,J)*(TG(I,J)-TP(I,J))*HTC 
END  IF 

M  ARTIFICIAL  VISCOSITY  ** 


ARVISR(I,J)=ARV*(ABS(AVG(UG,I+1,J,2,1)+0.1)+2.74E5) 
ARVISZa,J)=ARV*(ABS(AVG(WG,I,J+1,2,2)+0.1)+2.74E5) 
AHV  Ii.S(ENDR+1  ,J)=ARVISR(ENDR, J) 

ARVISR( I , ENDZ+ 1 )=ARVISR( I ,  ENDZ) 

ARVISZ(ENDR+1 , J)=ARVISZ(ENDR, J) 

ARVISZ{ I ,ENDZ+1 )=ARVISZ( I , ENDZ) 

”  PRIMARY  VARIABLE  SHIFT  M 

DO  30  K=3,2,-1 
RH01( I , J,K)=RH01 ( I , J,K- 1 ) 

RH02( I , J ,K)=RH02( I , J ,K- 1 ) 

UG( I , J , K) =UG( I , J , K- 1 ) 

WG(I,J,K)=WG(I,J,K-1) 

UP(I, J,K)=UP(I,J,K-1 ) 

WP(I,J,K)=WP(I,J,K-1) 

EG(I, J,K)=EG(I, J,K-1 ) 

EP(I,J,K)=EP(I,-J,K-1) 

CONTINUE 

CONTINUE 

CONTINUE 

RETURN 

END 


c 


SUBROUTINE  SWEEP 


“w***V\V*. 


C 

c 

c 


c 

c 

c 


c 

c 


c 

c 

c 

c 
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«  LOCAL  VARIABLES  M 

REAL  CGI,  CG?,  CG3 ,  CG4,  CPI,  CP2,  CP3,  CP4 

REAL  MCR1 ,  MCR2,  MGR3,  MGE4,  MGR5,  MGR6 ,  MGR7,  KGRtJ,  MGR9 

REAL  HGZ1 ,  MGZ2,  MCZ3,  HGZ4,  MGZ5,  MGZ6,  MGZ7 ,  HGZ8,  NGZ9 

REAL  MPR1 ,  KPR2 ,  MPR3,  MPR4,  MPR5,  MPR6 

REAL  MPZ1 ,  MPZ2,  MPZ3,  MPZ4,  MPZ5,  KPZ6 

REAL  EG1 ,  EG2,  EG3,  EG4,  EG5,  EG6,  EG7,  EG8,  EC9,  EGIO 

REAL  EP1 ,  EP2,  BP3,  EP4,  EP5,  EP6 

REAL  EGCHEM,  EPCHEM 

REAL  UGRHOI ,  WGRHOI ,  UPRH02,  VPRH02,  EGRHQ1,  EPRH02 
REAL  TAUGGR(0:6,0:6l),  TAUGGZ(0:6,0:6l ) 

INTEGER  I,  J 


**  COHHON  VARIABLES  ** 


REAL  RH01<0:6, 0:61,3),  RH02(0:6,0:61 ,3) ,  UG(0:6,0:6l ,3), 

A  UP{0:6,0:61,3),  WG(0:6,0:61,3),  WP(0:6,0:61,3), 

A  EG(0:6,0:6l ,3) ,  £P<0:6,0;6l , 3) 

REAL  RHOG(0:6,0:6l),  RHOP(O:6,0:6l),  PG(0:6,0:6l),  PP(0:6,0:6l), 

A  PGPHI(0:6,0:6l ) ,  PP1KF<0:6,0:61 ),  TG(0;6,0:6l),  TP(0:6,0:61 ), 

A  PHI(0:6,0:6l ) ,  RADP(0:6,0:6l ) ,  MDG{0:6,0:61), 

A  GAMMAC(0:6 ,0:61 > ,  DRAGR(0:6,0:61),  DRAGZ(0:6,0:61), 

A  QD0T(0:6 ,0:61 ) ,  ARVISR<0:6,0:61),  ARVISZC0:6,0:61) 

INTEGER  IGN(0:6,0:61) 

REAL  RLEN,  ZLEN,  DELR,  DELZ,  DELT,  STEP,  EPS,  ARV,  TOL,  R,  TIME, 

A  CVG,  CVP,  RO,  MW,  B1,  TGO,  TPO,  TCH,  TIGM,  B,  BUT,  EPT,  RHOPO, 

A  PHIO,  PGO,  PCH,  ?WALL,  RADPO,  KO,  HUGO,  PR 
INTEGER  INIT,  TOGGLE,  ENDR,  ENDZ,  STOPIT,  COUNT,  MAX, 

A  RPR I NT,  ZPRINT ,  TPRINT,  NTGR,  NTGZ,  MTPR,  NTPZ,  NPGR,  NPGZ 
COMMON  /PRI/  RH01,  RH02,  UG,  UP,  WG,  WP,  EG,  EP 
COMMON  /SEC/  RHOG,  RHOP,  PG,  P?,  PGPHI,  PP1MF,  TG,  TP,  PHI, 

A  RADP,  IGN,  MUG,  GAMKAC,  DRAGR,  DRAGZ,  QOOT,  ARVISR,  ARVISZ 
COMMON  //  RLEN,  ZLEN,  DELR,  DELZ,  DELT,  STEP,  EPS,  ARV,  TOL,  R, 

A  TIME,  CVG,  CVP,  RO,  MW,  B1,  TGO,  TPO,  TCH,  TIGN,  B,  BUT,  EPT, 

A  RHOPO,  PHIO,  PGO,  PCH,  PWALL,  RADPO,  KO,  HUGO,  PR,  INIT,  TOCCIE, 
A  ENDR,  ENDZ,  STOPIT,  COUNT,  MAX, 

A  RPRIHT,  ZPRINT,  TPRINT,  NTGR,  NTGZ,  MTPR,  NTPZ,  NFGR,  NPGZ 


EXTERNAL  AVG,  DBLAVG 
EGCHEMsQ . 95*EPT 
EPCHEM=0.05#EPT 
DO  10  1=1, ENDR 
DO  20  J= 1 , ENDZ 


*  GAS  CONTINUITY  EQUATION  • 

HHIftfHttHIHiHIHIIH 
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R=DELR*( 1-0.5) 

CG1=(AVG(RH01,I+1,J,2t1)*(R+D£LR/2.0) 

&  *UG(I+1 , J,2)-AVG(RH01 , I, J,2, 1 )*<R-D£LR/2.Q) 

&  »UC(I,J,2))/D£LR 

C 

CG2=(AVG{RH01 ,1 , J+1 ,2,2)*VGf I, J+1 ,2) 

&  -AVG(RH01 , 1,J, 2,2)eKG(i ,»)  2))/DELZ 

C 

CG3=(A»G(ARVISR;I+1,J, 1 , 1 )»(RB01( 1+1 , J,2)-RH01(I , J,2) ) 

&  -AVG'ARVISR, I , J, 1 , 1  )*(RH01{ I, J  ,2)-RB01{I-1 , J,2)))/DELR 

CG3=0.0 

C 

CGM=( AVG{ ARVISZ , I , J+ 1  1 ,2)*(RH01 ( I , J+1 ,2)-RH01 ( I , J,2) ) 

A  -AVGCARVISZ, I ,  J, 1 ,2)*(RH01 ( I , J ,2)-RH01 ( I , J-1 ,2)) )/DEL2 

IF(J.EQ.I)  CG4=(AVG(ARVISZ,I , J+2, 1 , 2)*(RH01(I , J+2,2) 

A  -RH01( I ,  J+1 ,2) )-AVG( ARVISZ, I , J+1 , 1 ,2)*(RH01(I, J+1 ,2) 

A  -RH01(I,J,2)))/DEL2 

C 

RH01 ( I , J, 1 )=RH01 (I, J,3)+2. 0*DELT*(-1 .0/R*CG1-CG2+CG3+CG4 
A  +GAMMAC(I, J)) 

C 

C  flHHHHMHIHlHtmiHfHtl 

C  *  PARTICLE  CONTINUITY  EQUATION  * 

C 

CP1=(AVG(RH02,I+1,J,2,1)*(R+DELR/2.0) 

A  *UP( 1+1 ,  J,2)-AVG(RH02, 1 , J,2, 1 )*(R-DELR/2.0) 

A  *UP(I,J,2))/DELR 

C 

CP2=(AVG(RH02,I,J+1,2,2)*VP(I,J+1,2) 

A  -AVG(RHQ2, I,J,2,2)*WP( I, J,2))/DELZ 

C 

CP3=(AVG<ARVISR,I+1,J,1,1)«(RH02(I+1,J,2)-RH02(I,J,2)) 

A  -AVG(ARVISR,  I,  J,  1 1  )*(RH02(  I ,  J  ,2)-RH02(I-1 ,  J,2) )  )/D£LR 

CP3=0.0 

C 

CP4=( AYG( ARVISZ, I , J+1 , 1 ,2)*(RH02( I, J+1 ,2)-RH02( I , J,2) ) 

A  -A7G(ARVISZ, I , J, 1 ,2)*(RH02( I , J ,2)-RH02( I , J-1 ,2) ) )/DELZ 

IF(J.EQ.I)  CP4=(AVG( ARVISZ, I ,J+2, 1 ,2)*(RHD2(I , J+2,2) 

A  -RH02C , J+ 1,2) )-AVG(AR?ISZ, I , J+1 , 1 ,2}«(RH02(I , J+t ,2) 

A  -RH02( I , J,2) })/DELZ 

C 

RH02 ( I , J , 1 ) = RH02 ( I , J , 3 ) +2 . 0*DELT*( - 1 . 0/R«CP 1 -CP2+CP3+CP4 
A  -GAHMACd.J)) 

C 

CJUJUUUUUI  ALA JUUI  ftAAtt  AJUUI MJIMJI 

RVVIRVVIVIVRVVVRffVVVVVVVV 

C  *  GAS  MOMENTUM  EQUATION  * 

C  * _ RADIAL  DIRECTION  • 

C 

R=DELR*( 1-1 ) 

IF(R.EQ.O.O)  GO  TO  1000 
HGR1=(RH01(I,Jf2)»(R+D£LR/2.0) 

A  •(AVG(UG,I+1,J,2,1))**2.0-RH01(I-1»J,2) 
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o  o  n  n  o  o 


4  »(R-DELR/2.0)*(AVG(UG,I,J,2,1))**2.o)/DELR 

«CR2s ( DBLAVG( RHO 1 , I , J+ 1 , 2 ) 

4  vAVG(UG, I , J+1 ,2t2)*AVG(WG,I , J+1 ,2,1) 

4  -DBLAVG{RH01,I,J,2)»AVG(UC,I.J,2,2) 

4  #AVG(WG, I ,J,2, 1 ) )/DELZ 

IF(TOGGLE.EQ.I)  THEN 

MGR3=(((R+DELR)*UG(I+1,J,2)-R*0G(I,J,2))/{R+DELS/2.0) 

4  -(R#UG{ I, J , 2)-(R-DELR)*UG(I-1 , J,2) )/(R-DELR/2,Q ) ) 

4  /(DELR**2.0) 

MGRMUG(I,J+1,2)+UG(I,J-1,2)-2.Q«UG(I,J,2))/D£LZ**2.0 

MGR5=NCR3+( WG( I ,  J+1 ,2)-VG( I , J,2)-(WG( 1-1 , J+1 ,2) 

4  -WG( 1-1 , J,2) ) )/(DELR#DEL2) 

ELSE 

MGR3=0.0 

MGR4=0.0 

HGR5=0.0 


END  IF 
C 

MGR6=(GAMMAC(I,J)+GAMMAC(I-1 ,J))/2.0*UP(I, J,2) 

MGR7  = (PGPHI ( I , J ) -PGPHI (1-1 »J)) /DELR 

MGR8=(ARVI5R(I, J)*(A¥G(R)851 ,1+1 , J,2, 1 )*UG(I+1 , J ,2) 
4  -AVG(RH01 , I , J,2, 1 )*UG{ I , > ,2) )-ARVISR( I~1 , J) 

4  *{AVG(RH01 , I, J,2, 1 )*UG( I , J,2)-AVG(RJfQ1 ,1-1 , J,2,2) 
4  *OG( I- 1 , J , 2 ) ) ) /DELR 

MGR9=(DBLAVG( ARVISZ, I, J+1 , 1 )*(AVG(RH01 , I , J+1 ,2, 1 ) 

4  *UG( I, J+1 ,2)  .,'VG{RH01 , 1 , J,2, 1  )*UG( I ,  J,2) ) 

4  -DBLAVG( ARVISZ , I , J, 1 )*( AVG(RH01 , I, J,2, 1 )*UG( I , J,2) 

4  -AVG(RH01,I,J-1,2,1)*UG(I,J-1,2)))/DELZ 

UG{I, J, 1 )=(AVG(RH01 ,I,J,3,1)*IK»(I,J,3) 

4  +2 ,0*DELT*( - 1 . 0/R»MGR 1 -HGR2+MUG( I , J ) «(KGR3+MGR4 

4  +1.0/3. 0*MGR5 ) +MGR6-DRAGR( I , J ) -MGR?  +MGR8+MGR9 ) ) 

4  /A¥G(RK01 ,1, J, 1,1) 


*  GAS  MOMENTUM  EQUATION  * 

*  AXIAL  DIRECTION  » 


1000  R=DELR*( 1-0.5) 

MGZ1s(DBLA¥G(RH01 , 1+1 , J,2)*(R+DELR/2.0) 
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'Cte  ^  —  ^*4  *.  V  mm9 


a. 


non 


A 

A 

A 


A 

C 

c 


*A¥G(aG,I+1,J,2,2)*AVG(VG,r+t,J,2tt) 
-DBLAVG(RH01 ,1 , J,2)*(R-DELH/2.0) 

#AVG(UG,  I ,  J,2,2)*AVG(t#G,  I ,  J  ,2, 1 ) )/D£LR 

-RHOKI.J- 1 ,2)*{AVG(WG,I ,J,2,2) )**2.0)/D£LZ 
IF(T0GGLE.EQ.  1 }  THEN 
IF(I.EQ.EMDR)  THEN 


**  NO-SLIP  BOUNDARY  CONDITION  ** 


& 

& 

C 

c 


c 

c 

c 


A 

A 

C 

c 

c 

c 

c 

c 

c 

c 

c 


A 

A 

A 

C 

c 


MGZ3= ( ( R+DELR/2 . 0 ) # ( - 1 . 0*WG( I ,  J ,  2 ) 
-1.0/3.0*VG(I-1,J,2))-(R-DELR/2.0)»(WG(I,J,2) 

-MG(I- 1 ,J,2) ) >/(R*DELR**2.0} 

ELSE 

MGZ3=( (R+DELR/2. 0)*(WG(I+1 , J,2)-WG(I>J*2) ) 
-(R-DELR/2.0)#(HG(I,J,2)-NG(I-1, J,2)))/(R*DELR**2.0) 

END  IF 

MGZilrCWGd.J+I^J+WCa.J-I^J-a.O^hJGd.J^JJ/DQ^^.O 

HGZ5=( (R+DELR/2. 0)*(UG(I+1 ,J,2)-UG( 1*1, J-1, 2)) 

- ( R-DELR/2.0 }•< UG(I , J , 2 >-UG(I , J- 1 , 2) ) ) /( R*DELR*DELZ ) 
+MCZ4 

ELSE 


MGZ 3-0.0 
HGZ4=0.0 


MGZ5=0.0 
END  IF 

MGZ6 = ( GAMMAC (I , J } +GAMHAC (I , J- 1 ) ) /2 . 0»«P (I , J , 2 ) 

MGZ7= ( PGPHI (I , J ) -PGPHI (I , J- 1 ) ) /DELZ 
IF(I.EQ.ENDR)  THEN 

HCZ8=(DBLAVG(ARVISR,I,J,1)»(AVG(HH01,I,J,2,2) 

*WGd,J,2)-AVG(RH01,I-1,J,2,2)*WG(I-1,J,2)) 

-DBLAWG(AR¥ISR,I-1,J,1)*(A¥G(RH0M-t,J,2,2) 

*WG( 1-1 , J ,2)-AVG(RH01 , I-2,Jt2,2)*VG(I-2, J,2)))/DELR 

ELSE 

MGZ8=(DBLA¥G(AR¥ISR, 1+1 , J» 1 )*(A7G(RH01 , 1+1 ,J,2(2) 
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ft*  ft*  ft*  fl*  R*  R*  fi*  11*  ft*  ft*  ft*  ft*  ft*  ft*  ft*  ft*  ft*  ft*  ft*  ft*  ft*  ft*  ft* 


»WG( 1+ 1 , J . 2 ) - A VG( RHO 1 , I , J , 2 , 2 >*WG( I , J , 2 ) ) 

-DBLAVGC ARVISR, I ,  J,  1 )*( AVG(RBQ1 ,1, J,2,2)*WG{ I , J,2) 
-A¥G(RH01 , 1-1 , J,2, 2}#VG(I-1 , J,2) ) )/OSLR 

END  I? 

IF(J.EQ.I)  THEN 

MGZ9=( ARVISZC I , J+1 )*(A¥G<RHOl ,1, J+2,2,2) 

*WG(I, J+2 ,2)-A¥G(RH01 , I , J+1 ,2,2)rfWG(I , J+1 ,2) ) 
-ARVIS2(  I , J)*(A¥G(RH01 ♦ I , J+1 ,2,2)*WG(I, J+1 ,2) 
-AVG(RH01 , I , J,2 ,2)*WG( I , J,2) ) )/DELZ 
ELSE 

MGZ9=(ARVISZ(I,J)*(AVG(RHOt,I,J+1,2,2)»WG(I,J+1,2) 
-AVG(RH01 . I ,J,2,2)*WG(I, J,2) )-ARVISZ(I, J-7 ) 
*(A¥G(RH01 , I , J,2,2)*WG(I , J ,2)-A¥G(RH01 , I, J- 1 ,2,2) 
*WG( I , J- 1 , 2 ) ) ) /DELZ 
END  IF 

WC( I ,  J,  1 )=( A¥G(RH01 , I, J ,3»2)*NG(I, J,3) 
+2.0*DELT*(-1.0/R*MGZ1-MGZ2+HUG(I,J>»(MC£3+MCZ1» 

+ 1 . 0/3. 0*MGZ5 ) +MCZ6-DRAGZ (I, J) -HGZ7+MGZ8+HGZ9  > ) 
/AVG(RH01 , I, J, 1 ,2) 

*  PARTICLE  MOJffiNTW  EQUATION  * 

»  RADIAL  DIRECTION  * 

fHHHHHtHMHIMHIlHWtH 


R=DELR*( 1-1 ) 

IF(R.EQ.O.O)  GO  TO  2000 

MPR1=(RB02(I,J,2)»(R+DELR/2.0) 

*(A¥G(UP,I+1 , J,2, 1 ) )**2.0-RH02( 1-1 ,J, 2) 
*(R-DELR/2.0)*(AVG(UP, I# J ,2, 1 })**2.0)/DELR 

MPR2=(DBLAVG(RH02, I , J+1 ,2) 
*A¥G(UP,IfJ+1,2,2)*A¥G(HP,I,J+1 ,2, 1) 

-DBLAVG(RH02, I , J,2)*A¥G(UP, I , J,2,2) 

*A¥G(VP,I,J,2, 1))/DELZ 

MPR3= ( GAMMAC ( I , J ) +GAHMAC ( I - 1 , J ) ) /2 . 0*UP( I , J „ 2 ) 

MPR4=(PP1MF(I,J)-PP1MF(I~1 ,J))/DELR 

MPR5=(AR¥ISP(I,J)*{A¥G(RJS)2,I+1,J,2, 1  )*UP( 1+1 , J, 2) 
-AVG(FH02, I , J,2, 1 )*UP( I , J ,2) )-ARVISR( 1-1 ,J) 
*(A¥G(RH02,I,J,2,1)*UP(I,J,2)-A¥G(RH02,I-1,J,2,1) 
*UP( 1-1 , J,2) ) )/DELR 

MPR6s(DBLA?C<ARVISZ, I, J+1, 1)»(AVG(RH02,I,J+1, 2,1) 
*UP( I, J+1 ,2)-AVG{RH02, I,J ,2,1 )#OP(I , J,2)) 
-DBLAVG(ARVISZ,I,J,1)*(A¥G(KH02,I,J,2f1)*UP(I,J,2) 
-A¥G(RH02, 1 ,  J- 1 ,2, 1  )*UP(  I  ,J- 1 ,2)  ))/DOZ 

UP(I,J,1)=(A¥G(RB02,I,Jf3,1)*OT(I.J,3) 
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4 

4 


+2 .Q*DELT*(- 1 ,0/R*MPR1-MPR2“MPR3+DRACR(  I ,  JJ-MPR4+HPR5 
+MPR6))/AVG(RH02,I,J,1,1> 


C 

C 

C 

C 

C 

C 

2000 


4 

4 

4 

C 

4 

C 

C 

C 

& 

4 

4 

C 


4 

4 

4 


4 

4 

4 


•  PARTICLE  MOMENTUM  EQUATION  * 

*  AXIAL  DIRECTION  * 


R=DELR*(I-0.tU 

MPZ 1  =  ( DBUl VG ( RH02 , 1 ♦ 1 , J , 2 ) * ( R+DELR/2 . 0 ) 

*AVG(UF,  I-*- 1 , J,2,2 )*AVG(NP, 1+1 ,  J ,2, 1 ) 

- DEL AVG( RH02 , I , J , 2 )»( R-DELR/2 . 0 > 

*AVG(UP, I , J,2,2)*A¥G(WP, I ,J,2, 1 )}/DELR 

MPZ2=(RH02( I , J,2)*( AVG{WP,I , J+1 ,2,2) )**2.0 
-RH02( I , J-1 ,2)*( AVG(VP, I , J,2,2) )**2.0)/DELZ 

MPZ3=  ( GAMMAC  ( I ,  J  )+GMHAC  ( I- 1 ,  J- 1 ) )  /2 . 0*WP(  I ,  J ,  2 ) 

MPZ4=(PP1MF( I , J)-PP1MF( I ,J-1 ) )/D£LZ 

MPZ5={DBLAVG( ARVISR, 1+1 , J, 1 )#(A¥G(RHQ2, 1+1 , J,2,2) 
»MP(I+1,J,2)-AVG(RH02,I,J,2,2)«VP(I,J,2)) 

-DBLA?G( AR? ISR , I , J , 1 )*( AVG( RH02 , I , J , 2 , 2 )«WP( I , J , 2 ) 
-AVG{RH02, 1-1 , J,2,2)*NP(I-1 ,J,2)))/D£LR 

IF(J.EQ.I)  THEN 

MPZ6=(ARVISZ{ I, J+1 )*{A¥G(RH02,I , J+2,2,2) 

*VP(I , J+2,2)-AVG(RHD2,I , J+1 ,2,2)*WP(I, J+1 ,2)) 
-AR¥ISZ( I , J)*(AVG(RH02,I , J+1 t2,2)*WP(I,J+1 ,2) 

-  A  VG  ( RH02 , 1 ,  J ,  2 , 2 )  *W  ( I ,  J ,  2 ) ) )  /  DELZ 
ELSE 

HPZ6=(AR¥ISZ(I,J)«(AVG(RH02,I,J+1,2,2)«WP(I,J+1,2) 
-AVG(RH02,I , J,2,2)*HP(I, J,2))-AR¥ISZ(I,J-t } 
•(AWGlRHOZ.I.J.Z^JWd.J^J-AyGtR^.I.J-l^^) 
•VP(I,J-1,2)))/DeLZ 
END  IF 


C 

WP(  I ,  J,  1  )=(AVG(RH02, 1 ,  J  ,3.2)*W(I,J  ,3) 

4  +2.0*DELT*{ - 1 . 0/R*MPZ 1 -MPZ2-MPZ3+DRACZ( I , J )-MPZ4+HPZ5 

4  +MPZ6) )/A?G(RH02, I , J, 1 ,2) 

C  mm mm,, 

C  «  GAS  ENERGY  EQUATION  • 

c 

EG1=( AVG(RH01 , 1+1 ,J,2, 1 )*UG(I+1 , J,2)*(R+DELR/2.0) 

4  *(A¥G(EG,I+1,J,2,1)+((PGPHI(I+1,J)+PCPHI(I,J)}/2.0) 

4  /AVG(RH01,I+1,J,2,1))-AVG(RHD1,I,J,2,1) 

4  *OG(I , J,2)*{R-DELR/2.Q)*(A¥G(EG,If J,2, 1 ) 

4  +((PGPHI(I,J)+PGPHI(I-1,J))/2.0) 

4  /A?G{RBQ1,I,J,2,1)))/DELR 

C 

EG2=(A¥G(RHD1,I,J+1,2,2)»NG(I,J+1,2) 
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& 

& 

& 

& 

& 

c 

c 

c 

c 

c 

& 

& 

& 

c 

c 

c 

c 


•( AVG(EG, I , J+1 , 2t2)+(  (PCPHKl ,  J+t  )+PCPHI(  I,  J))/2.0) 
/AVC(R!©1 ,1,  J+l  ,2.2) )-AVC(RHQ1 ,1 ,  J,2,2) 
*WCd,J.2>*(M?G(EG,I,.J.2.2> 
+((PGPHI(I,J)+PCPHI(I,J-1))/2.0) 
/A¥G{RH01fI,J,2,2)))/DELZ 

IF(TOGGLE.EQ.I)  THEN 

EG3= (UG( 1+ 1 , J , 2 )-UC( I , J , 2) ) /DEL* 

EG4=AVG(UG, 1+1 , J ,2, 1 )/R 

EG5= ( WG( I , J+ 1 , 2 )- WG( I , J ,2 ) ) /DELZ 

EG6= < DBLAVG( NG , 1+1 , J+ 1 , 2 ) 

-DBLAVG(WG,I , J+1 ,2) )/DELR 
+(D8LAVG(UG, 1+1 , J+1 ,2) 

-DBLA ¥G ( UG , I + 1 , J , 2 ) ) /DELZ 

ELSE 

EG3=0.0 

EG4=0.0 

EG5=0.0 


EG6=0.0 


END  IF 

EG7=  ( DRAGR(I+ 1 ,  J  )+DRAGR(I ,  J )  )•(  UP(  I*  1 ,  J , 2  )+tJP(  I ,  J ,  2 ) )  /*  .0 
+(DRAGZ(I,J+1 )+DRAGZ(I tJ))*(UP(I,J+1,2)+ilP(ItJl2))/4.0 

EG8=GAMMAC( I , J)*( EGCHEM+A¥G(UP, 1+1 , J,2, 1 )**2. 0/2.0 
+AVG(WP, I , J+1 ,2,2 )**2. 0/2.0) 

EG9=(AVG(AR?ISR, I*' fid . ’ )*(RH01{ 1+1 , J,?)*®!! 1+1 . J,2) 

-RH01 ( I , J , 2 )*EG( I , J , 2 ) )-AVG( ARVISR, I , J , 1 , 1 )»{RH01( I ,J , 2 ) 

•EG ( I , J , 2 ) - RH0 1 ( I - 1 , J , 2 } •EG ( I - 1 , J , 2 ) ) ) / OELR 
IF( I.EQ. 1 } 

EG9=(A?G(ARVISR,I+2, J, 1 , 1)*(RH01{I+2,J,2)*EG(I+2,J,2) 

-RH01(I+1,J,2)«EG(I+1,J,2})-A?G(ARVISRtI+1tJd»5) 

•(RH01( 1+1 , J*2)#EG( 1+1 , J»2)-RH01(I,J»2)*BG(I,Jt2}))/DELR 
IF(I.EQ.ENDR) 

EG9=(A?G(ARVISR,I,*Id.D*(RH01(I,J*2)*K(I,J,2) 

-RH01( 1-1 , J*2)*EG( 1-1 ,J ,2) )-AVG(ARVISR, 1-1 . J ,1,1) 
•<RHoid-i,j,2>*EGd-id,2)-w»iu-2,j,2)*8cd-2,J.2)))/DiLR 

EG10=<A¥G(AR¥ISZ,  I ,  J+1 ,1  t2)»<RHD1d ,  J+1 ,2>«EG(  I ,  J+1 ,2) 

-RHOKI ,  J»2)*EGd ,  J  ,2)  )-A¥G(  ARVISZ,  I  ,«!*  1  *2)*(RBQ1d ,  Jf2) 
*EG(I,J,2)-RH01d,J-1,2)*EG(I,J-1,2)))/D€LZ 
IF(J.EQ.I)  K!1G=( AVG(ARVISZ, I » J+2, 1 ,2)*(RH01( I, J+2»2) 

•EG(I, J+2,2)-RHD1(I,J+5  »2)#EG(I,J+1 12)) 
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C 

C 

c 

c 

c 


c 

c 


l 

t 


t 

& 

& 

& 

& 


& 

& 

& 

& 

& 


& 

& 

& 

& 

A 


& 

& 


& 

& 

& 


& 

& 

& 


& 

& 


A 

& 

& 


-AYG( ARVISZ, I » J+1 » !,2)*(RHD1(I, J+t ,2)*£G(I,J+1,2) 
-RH01(I,J,2)#EG(I,J,2)))/DELZ 


EG(I,J,1)=(RH01(I,J,3)*EG(I,J,3)+2.0*DELT*(-1.0/R*EG1-EG2 
+HUG( I , J )*( 2 . 0* ( EG3**2 . 0+EG4**2 . 0+EG5**2 . 0 ) 

-2 . 0/3 . 0*(  EG3+EG4+EG5 )**2 . 0+EG6**2 . 0+A  ?G(  UG ,  1+ 1 ,  J ,  2 , 1 ) 

*  ( MGR3+MGR4+ 1 . 0/ 3 . 0*MGR5 )  +A ?G  ( WG ,  I ,  J+ 1 , 2 , 2 ) 

* ( MGZ  3+MGZ4+ 1.0/3. 0*MGZ5 ) ) 

-QD07 (I, J) — EG7 +EG8+EG9+EG 10)) /RHO 1  ( I ,  J ,  1 } 


•  PARTICLE  ENERGY  EQUATION  * 

*•***»*«*#*«»**#**«***»««•** 


EP1=( AVG(RH02, 1+1 , J ,2 , 1 )*UP( 1+1 , J ,2)*(R+DELR/2.0) 
*(AWG(EP,I+1 , J,2, 1 )+( (PP1MF(I+1 t J)+PP1MF(I, J) )/2.0) 
/AVG(RK02, 1+1 , J,2, 1 ) )-AVG(RH02, 1 , J,2, 1 ) 

*UP{ I , J,2)*(R-DELR/2.0)*(AVG(EP, I, J,2, 1 ) 

+( (PpiMPd ,  J)+PP1MF(I-1 ,  J )  )/2.0) 
/AVG(RH02,I»J>2,1)))/DELR 


EP2=(AVG(RH02,I,J+-1»2,2)«WP(I<J+1.2) 

*(AVG{EP, I • J+1 *2,2)+( (PPIMf ( I , J+1 )+PP1MF(I, J) )/2.0) 
/AVG(RH02r I i J+1 , 2>2) )-AVG(RH02, I • J»2,2) 

*WP{ I ,J»2)*(AVG(EP, I , J.2,2) 
+({PP1MF(I,J)+PP1MF(I,J-D)/2.0) 

/A VG { RH02 , I , J , 2 , 2 ) ) ) /DELZ 


EP3=EG7 


EP4=GAMMAC( I , J )*( EPCHEM- AVG( UP , I + 1 , J , 2 , 1 ) »*2 . 0/2 . 0 
-AVG<VP,J;,J+1,2,2)«*2. 0/2.0/ 


EP5=(AVG(ARVISR,I+1|J,1,D*(RH02(I+1,J,2)«EP(I+1,J,2) 

-RH02(I, J,2)*EP( I , Jf2) )-AVC(ARVISR, I, J ,1*1 )*(RHQ2( I, J,2) 

*EP( I ,J»2)-RH02(I-1 , J,2)*EP(I-1 , J,2)))/DELR 
IF( I , EQ. 1 ) 

EP5=(AVG(ARVISR,I+2,J,1,1)*(RH02(I+2,J,2)»EP(I+2,J,2) 
-RH02(I+1 , J*2)*EP( 1+1 , J)2) )-A?G( ARVISR, 1+1 1  *  1 ) 

*(RH02( 1+1 , J,2)*EP(I+1 , J|2)-RH02(I , J,2)*EP(I , J,2) ) )/DELR 
IF( I .EQ.ENDR) 

EP5=(AVG(ARVISR,I,J,1,1)*(RH02(I,J,2)*EP(ItJ»2) 

-RH02( 1-1 , J,2)*EP( I- 1 , J  »2 ) )-AVC( ARtf ISR, I— 

#(RHG2( 1-1 , J , 2 )*EP ( I— 1 , J»2)-RH02(I-2, J  »2)*EP(I-2, J,2) ))/DELR 


EP6=(AVG(AP.VISZ,I,J+1,1,2)*(RH02(ItJ+1,2)*EP(I,J+1,2) 

-RH02( I , J? 2)*EP( I , J  »2) )-AVG( ARVISZ, I , J 1 1 *2)*(RH02( I , J,2) 
*EP( I , J , 2 )-RH02( I , J- 1 ,2 )»EP( I , J- 1 , 2 ) ) ) /DELZ 
IF(J.EQ.I)  EP6= ( AVG ( ARV ISZ , I , J+2 , 1 , 2 )*( RH02( I , J+2 , 2 ) 

*EP( I , J+2,2)-RH02( I ,  J+1 ,2)*EP(I , J+1 ,2) ) 
-A?G(ARVISZ,I,J+1,1,2)*(RH02(I,J+1,2)*EF(I,J+1,2) 
-RH02(I, J»2)*EP(I , J*2) ) )/DELZ 


EP(I , J , 1 )s(RH02( I, J»3)*EP( I i J,3)+2.0*DELT*(-1 .0/R*EP1-EP2 
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c 

c 


„c£ 
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200 

>> 
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300 

400 
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600 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 


+QDOT { I , J)  +EP3+EP4+EP5+EP6 ) ) /RH02 (I , J, 1 ) 


IF(COUNT.GE. 00. AND. COUNT. LE. 10}  THEN 
WRITE(6 , 100)  I,J 

WRITE(6 ,200)  CG1,CG2,CG3fCG4,(RH0l!l,J,K),K=1,3) 
WRITE(6,200)  CP1,CP2,CP3,CP4, (RH02(I,J,K) ,K=1 ,3) 

WRITE (6, 300)  HGR1,MCR2,HUG(I,J),MGR3,HGR4,MGR5,MGR6, 

4  DRAGR(I,J),MGR7,MGR8,MGR9,(UC(I,J,K),Kr1,3) 

WRITE(6,300)  MGZ1,MGZ2,MUG(I,J),MGZ3,MGZ4,MGZ5,MGZ6, 

&  DRAGZ(I,J),MGZ7,MGZ8,MGZ9,(WG(I,J,K),K=1,3) 

WRITE(6 ,400 )  MPR1,MPR2,MPR3,HPR4,MPR5,(UP(I,J,K),K=1,3) 
WRITE! 6, 400)  MPZ1,MPZ2,MPZ3,MPZ4,MPZ5,(WP(I,J,K),K=1,3) 
WRITE(6 ,500)  EG1,EG2,EG3,EG4,EG5,EG6,EG7,EG8,EG9,EG10, 

A  (EG(I,J,K),K=1,3) 

WRITE(6,600)  EP 1 , EP2 , EP3 , EP4 , EP5 , EP6 , ( EP( I , J , K ) , K= 1 , 3 ) 
FORMAT!/, 'POSITION  ', IX, 2(13)) 

FORMAT(1X,7(2X,G10.4>) 

FORMAT{ 1X,10(2X,G10.4)/1X,4(2X,G10.4)) 

FORMAT! 1X,8(2X,G1Q.4)) 

FORMAT! 1X,10(2X,G10.4)/1X,3(2X,G10.4)) 

FORMAT! 1X,9(2X,G10.4)//) 

END  IF 

TAUGGR ( I , J ) sMUG! I , J ) *(MGR3+MGR4+ 1.0/3- 0*MGR5 ) 

TAUGGZ ( I , J ) =  MUG ( I , J ) * ( MGZ  3+MGZ4+ 1.0/3. 0»MGZ5 ) 

20  CONTINUE 

10  CONTINUE 

CALL  DUMP2 

*  BOUNDARY  CONDITIONS  * 

ffllfffllrlWffffff  ififf* 

«  AT  Z=0  END  WALL  SYMMETRY  CONDITIONS  ** 

DO  30  I=0,ENDR+1 
RH01( I ,0, 1 )=RH01(I, 1,1) 

RH02( 1,0,1 )=RH02( I ,1,1) 

UG( 1,0,1 )=UG( 1,1,1) 

WG!I,1,1)=0.0 
WG( 1,0,1 )=WG( 1,2,1) 

UP! i,op i )=up(i , 1,1) 

WP!I,1,1)=0.0 

WP( I ,0, 1 )=WP( 1,2,1) 

EG! 1,0,1 }=EG( 1,1,1) 

EP(I,0, 1 )=EP( 1,1,1) 

30  CONTINUE 

**  AT  R=0  CENTER  LINE  SYMMETRY  CONDITIONS  M 
DO  40  J=1,ENDZ 
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RH0I{O,J, 1)=RH01( 1,J,  1) 
HH02(0, J , 1 )=RH02( 1 3Jf 1 ) 
UG(1,J,1)=0.0 
UG(0, J, 1 )=UG(2,J, 1 ) 
WG(0, J, 1 )=WG( 1 ,  J,  1 ) 
UP(1,J,1)=0.0 
UP(0,J, 1 ):UP(2,J, 1 ) 
WP(0, J, 1 )=WP( 1 , J, 1 ) 
EG(0,J,1)=EG(1,J,1) 
EP(0, J, 1 )=EP( 1 , J, 1 ) 
CONTINUE 


**  AT  Z=ENDZ  END  WALL  RADIATIVE  CONDITIONS  « 


DO  50  I=0,ENDR+1 

RHO 1 ( I , ENDZ + 1 , 1 ) : { 5 • Q*RH0 1 ( I , ENDZ , 1 ) -4 . O*RH0 1( I , ENDZ- 1 , t ) 
k  4-RH0 1(1,  ENDZ-2 , 1 ) )  /2 . 0 

RH02 ( I , ENDZ* 1 , 1 ) = ( 5 . 0*RH02 ( I , ENDZ , 1 )-4 . 0*RH02( I , ENDZ- 1 , 1 ) 
k  +RH02(I,ENDZ-2, 1) )/2.0 

DG( I .ENDZ+1 , 1 )=(5.0*UG( I ,END2, 1 )-4.0*UG( I, ENDZ- 1 , 1 ) 
k  +0G( I, ENDZ-2, 1)}/2.0 

»G( I , ENDZ+ 1 , 1 ) = ( 5 . 0»WG ( I , ENDZ , 1 ) -4 . 0*WG(I , ENDZ- t , 1 ) 
i  +WG(  I ,  ENDZ-2 , 1 ) )  /2 . 0 

UP(  I , ENDZ-t- 1 , 1 )  =  ( 5 . 0*UP ( I , ENDZ ,  1 ) -4.0*UP( I , ENDZ- 1 , 1 ) 
t  +UP( I , ENDZ-2, 1 ))/2.0 

WP( I .EHDZ+1 , 1 )=(5*0#NP(I,ENDZ, 1 )-4.0*WP(I,ENDZ-l , 1) 
t  +WP{  I, ENDZ-2, 1 ))/2.0 

EG(I,ENDZ+1,1)=<5.0*EG(I,ENDZ,1)-4.0»EG(I,ENDZ-1,1) 
t  fEG( I, ENDZ-2, 1))/2.0 

EP(I .ENDZ+ 1 , 1 )= (5.0«EP( I , ENDZ, 1 )-4.0»EP(I,ENDZ- 1 , 1 ) 
i  +EP( I, ENDZ-2, 1))/2.0 

CONTINUE 


AT  RcENDR 


CIRCUMFERENTIAL  NALL  NO-SLIF  CONDITIONS  » 
CRACK  OPENING  CONDITIONS  t* 


NOTE:  NO-SLIP  CONDITIONS  ARE  APPLIED  ONLY  IF 

THE  VISCOUS  GAS  EFFECTS  ARE  TO  BE  INCLUDED. 
CONSEQUENTLY,  ACTUAL  BOUNDARY  CONDITIONS 
ARE  APPLIED  VIA  ONE  SIDED  DIFFERENCING  IN 
THE  MAIN  INTEGRATION  LOOP  IF  TOGGLE: 1. 

CONDITIONS  LISTED  HERE  ARE  ONLY  USED  TO  KEEP 
FALSE  EXTERIOR  NODES  FILLED  AND  HENCE  THE  VALUES 
MAY  OR  MAY  NOT  BE  CORRECT. 


DO  60  J= 1 , ENDZ 

RH01 ( ENDR+ 1 , J , 1 ) =RH01 ( ENDR , J , 1 ) 
RH02  ( ENDR-t- 1 ,  J ,  1 )  *  RH02  { ENDR ,  J ,  1 ) 
UG(ENDR+1 . J, 1 )=0.0 
WG{ENDR+1 ,J, 1 )=UG(ENDR,J, 1) 
UP(ENDR+1 , J, 1 )=0.0 
WP(ENDR+1 , J, 1 )=WG(ENDR, J, 1 ) 
EG(ENDR+1 , J, 1 )=£G(ENDR, J, 1 ) 
EP(EMDR+1 , J, 1 )sEP(ENDR, J, t ) 


— * tern •*4w**i 1 — 1 — - - ' — ^ 


1  . 


CONTINUE 


60 
C 
C 

C  *  REDUCTION  OF  TRUNCATION  ERRORS  * 

C  . . . . . 

c 

DO  70  I=0,£NDR+1 
DO  75  J=0,ENDZ+1 

IF( AB5(UG(I , J, 1)).LE.T0L)  UG(I,J, 1 )=0.0  A 

IF(ABS<VG(I,J,1)). LE.TOL)  NG(I, J, 1 )=0.0 
IF( ABS(UP( I ,  J, 1)) .LE.TOL)  UP(I, J, 1 )=0.0 
IF(  ABS( WP( I ,  J,  1 } ) .LE.TOL)  ffP<  I ,  J,  1  )=0.0 
75  CONTINUE 

70  CONTINUE 

C 

C  iinmiwmim 

C  *  TIME  SMOOTHER  * 

C  HVWHHMHHHHHHNI 

c 

DO  80  I=0,EMDR+1 
DO  85  J=0,ENDZ+1 

ECUUK)  1 = EG ( I , J , 2 ) • RH0 1 ( I , J , 2 ) +EPS* ( EG ( I , J , 1 ) *RH0 1 ( I , J f 1 ) 

&  -2.0*EG(ITJ,2)*RH01(I,Jf2)+EG(I,J,3)*RH01(I*J»3)) 

EPRH02=EP( I , J,2)»RH02( I , J ,2)+EPS»(EP( I , J , 1)*!W02(I , J, 1 ) 

&  -2.0*EP(I,J,2)»RH02(I,J,2)^EP(I,J,3)*RH02(I,J,3)) 

UGRH01=UG( I , J, 2)*RH01 ( I , J ,2)+EPS*fUG(I ,J, 1)*RH01( I, J, 1 ) 

4  -2.0*UG(I,J,2)»RH01(I,J,2)+UC<I,J,3>*RHO1(I,J,3>) 

NGRH01=NG( I , J,2)»RH01 ( I , J,2)+EPS«(NG( I , J, 1 )*RK01( I , J, 1 ) 

&  -2.0*WG{ I, J,2)#RH0t(I,J,2)+HG(I,J,3)*RHO1( I , J,3) ) 

UPRH02=UP< I , J , 2 )»RB02( I , J , 2 ) +EPS*( UP( I , J , 1 >*RHQ2( I , J , 1 > 

&  -2.0*UP(I,J,2)*RH02(I,Jt2)«4n*(I,J,3)*RH02(IfJ,3>) 

WPRH02=VfP( I , J ,2)*RH02( I , J ,2)+EPS*(NP( I , J , 1 )*RH02( I , J , 1 ) 

&  -2.0*HP(I, J,2)*RH02(I, J,2)+HP(I,J,3)*RB02(I,Jt3)) 

RBOKI, J,2)=RH01(If J,2)+EP5*(HH01(I, J, 1>-2.0*fiHD1{I,J,2) 

4  +RH01(IIJ,3)) 

RHQ2(I , J,2)=RHG2(I, J,2)+EPS*(RHQ2( I , J, 1 )-2.0*RH02(I , Jt2) 

4  +RH02( I , J ,3) ) 

UG{ I , J , 2 ) =UGRH0 1 /RB01 ( I , J , 2 ) 

VC( X, J,2)=WGRH01/RH01 ( I , J,2) 

UP ( I , J , 2 ) = UPRM02 / RH02 ( I , J , 2 ) 

HP(I,J,2)=HPRH02/RH02(X,J,2)  • 

EG( I , J , 2 ) sEGRBQI /RM01 ( I . J , 2 ) 

EP( I , J,2)=EPRH02/RH02( I ,J,2) 

IF<COUNT.GE. 00. AND. COUNT. LE. 10)  THEN 
WRITE(6,700)  I , J, RH01 ( I , J ,2 ) ,RH02{ I , J,2) ,UG{ I , J,2 ) ,  - 

4  WG( I , J ,2) ,UP(I,J,2),NP(I,J ,2) ,EG{ I, J,2) ,EP( I , J ,2) 

700  F0RMAT( IX, 2(13), 8(G 10. ft)) 

END  IP 

85  CONTINUE 

80  CONTINUE 
RETURN 
END 
C 
C 
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V> Vk V>  V>  Oi>  Vi  .I’AlUV.-.l,. CVAYVWVV'VUl.tS 


%. »  yN  17W  »  . ' "V ,v"fc  «  Tk.AJ 


non  o  o  o  ono  on  non  non  on  ooo 


C 

FUNCTION  AVG  (VALU,  SITER,  SITEZ,  TIME*  DIR) 

**  LOCAL  VARIABLES  M 

REAL  VALU(0:6, 0:61,3) 

INTEGER  SITER,  SITEZ,  TIME,  DIR 


IF  (DIR.EQ.1)  THEN 

AVG=(VALU(SITER, SITEZ, TIME)+VALU(SITER-1, SITEZ, TIME) )/2.0 
ELSE  IF  (DIR.EQ.2)  THEN 

AVG={VALU( SITER, SITEZ, TIME )+VALU(SITER, SITEZ- 1 .TIME) )/2.0 
END  IF 
RETURN 
END 


FUNCTION  DBLAVG  (VALU,  SITER,  SITEZ,  TIME) 

«  LOCAL  VARIABLES  ** 

REAL  VALU(0:6,0:6l,3) 

INTEGER  SITER,  SITEZ,  TIME 


DBLA VG= { V ALU<  SITER , SITEZ , TIME ) +VALU( SITER- 1 , SITEZ ,  TIIffi  ) 

&  +V ALU(SITER, SITEZ- 1 , TIME )+VALU(SITER-1 .SITEZ- 1 .TIME) )/4.0 
RETURN 
END 


SUBROUTINE  DUMP 

**  LOCAL  VARIABLES  *« 

REAL  RLOC,  ZLOC,  PGSI,  RADPSI ,  TIMOUT,  DTOUT,  PWALOT 

INTEGER  I,  J 

«  COMMON  VARIABLES  H 

REAL  RH01(Q:6,0:61,3),  RH02(0:6,0:6l,3),  UG(0:6,0:6l,3), 

A  DF(0:6, 0:61,3),  «G(0:6,0:6l ,3),  VP(0:6,0:6l ,3), 
k  EG(0:6,0:61,3).  EP(0:6,0:6l ,3) 

REAL  RHOG(0:6,0:6l),  RH0P(0:6,0:6l ) ,  PG(0:6,0:6l),  PP(0:6,0:61 ) , 
k  PGPHI(0:6,0:6l),  PP1MF(0:6,0:61) ,  TG(0:6,0:6l),  TP(0:6,0:61), 
k  PHI(0:6,0:61),  RADP(0:6,0:6l),  MUG{0:6,0:61 ), 
k  GAM<AC(0:6,0:61),  DRAGR(0:6,0:6l ) ,  DRAGZ(0:6,0:61), 
k  QDOT(0:6,O:61 ) .  ARVISR(0:6,0:61),  AR¥ISZ(0:6,0:6l) 

INTEGER  IGN(0:6,0:61) 

REAL  RLEN,  ZLEN,  DELR,  DELZ,  DELT,  STEP,  EPS,  ARf ,  TOL,  R,  TIME, 
&  CVG,  CVP,  RO,  MU,  B1,  TOO,  TPO,  TCH,  TIGN,  B,  BNT,  EPT,  RBOPO, 


103 


4  PHIO,  PGO,  PCH,  PWALL,  RADPO,  KO,  MUGG,  PR 
INTEGER  INIT,  TOGGLE,  ENDR,  EMDZ,  STOPIT,  COUNT,  MAX, 

&  RPR I NT,  ZPRINT,  TPRINT,  NTGR,  NTGZ ,  NTPR,  NTPZ,  NPGR,  NPGZ 
COMMON  /PRI/  RHOI ,  RH02,  UG,  UP,  MG,  VP,  EG,  EP 
COMMON  /SEC/  RHOG,  RHOP,  PG,  PP,  PGPHI,  PPtHF,  TG,  TP,  PHI, 

&  RADP,  IGN,  MUG,  GAMMAC,  ORAGR,  DRACZ,  QDQT,  ARVISR,  ARVISZ 
COMMON  //  RLEN,  ZLEN,  DELR,  DELZ,  DELT,  STEP,  EPS,  ARV,  TOL,  R, 

&  TIME,  CVG,  CVP,  RO,  HH,  B1,  TGO,  TPO,  TCHt  TIGN,  B,  BNT,  EPT, 

&  RHOPO,  PHIO,  PGO,  PCH,  PHALL,  RADPO,  KO,  HUGO,  PR,  INIT,  TOGGLE,  A 

&  ENDR,  ENDZ,  STOPIT,  COUNT,  MAX, 

&  RPRINT,  ZPRINT,  TPRINT,  NTGR,  NTGZ,  NTPR,  NTPZ,  NPGR,  NPGZ 
C 
C 

EXTERNAL  AVG 
PRINT  100 
DO  10  1=1, ENDR 

RLOC= I*DELR-DELR/2 . 0 
DO  20  J=1 ,ENDZ 
ZL0C=J*DELZ-DELZ/2 .0 
PGSI=PG( I,J)*1.0E-10 
RADPSI=RADP(I,J)*1.0E4 
UGPRT=AVG(UG, I, J,2, 1 ) 

WGPRT=AVG(WG, I , J,2,2) 

UPPRT=AVG(UP, I , J,2, 1 ) 

HPPRT=  AVG(  HP ,  I ,  J ,  2 , 2 ) 

IF  (( (J.EQ.O) .0R.(J.EQ.ENDZ).0R.( (J/ZPRINT)*ZPRINT.EQ. J)) 

4  .AND.((I.EQ.O).OR.(I.EQ.ENDR).OR.((I/RPRINT)*RPRINT.EQ.I))) 

&  PRINT  200 , RLOC , ZLOC , UGPRT , UPPRT , MGPRT , HPPRT , 

&  PGSI,TG(I,J),TP(I,J),RHOG(I,J),RHOP(I,J),PHI(I,J),RADPSI, 

&  GAMHACC I,J) ,QDOT(I,J},DRAGR(I,J) ,DRAGZ(I,J) 

20  CONTINUE 

PRINT  400 
10  CONTINUE 

PWAL0T=PWALL*1 .0E-10 
TIM0UT=TIHE*1.QE6 
DTOUT=DELT*  1 .0E6 

PRINT  300, TIMOUT, COUNT, DTOUT, PH ALOT 
C 
C 

100  FORMATCl'.lX.'RLOC'.SX.'ZLOC'.SX, ’UG',6X,'UP»,6X,'HG»,6X,'WP', 

&  3X,'GAS  PR£SS,,3X,'TG',5X,'TP*,4X,'RH0G*,3X,,RHDP,,3X,,PHI,f  « 

4  2X, 'RADIUS’ ,2X, ' GAPBfA * ,5X, 'QDOT' ,5X, 'DRAGR' ,3X, 'DRACZ'/ 

4  3X,2(’CH',5X),4( 'CM/S', 4X),'GPA',5X,2(’DEG  V ,3X), 'G/CC' , 

&  'n/rr'  in*  'fin'  /vkiil •  i 

200  FORMAT*'  ’ *2<FS.3, IX) ,4(F7.0, 1X),F8.5, 1X,2(F6.0, IX) ,2(F6.4, IX) ,  * 

4  F5.3,1X,F6.2,4(1X,G8.2)) 

300  FORMAT (/'  DATA  IS  RECORDED  AT  ',F8.4,'  MICRO-SECONDS' ,6X, 

4  15,'  INTEGRATIONS  AT  THIS  TIKE'//'  CURRENT  TIME  STEP  IS  *, 

4  F8.4, ’  MICRO-SECONDS’, 6X, 'MAXIMUM  HALL  PRESSURE  IS  ', 

4  F8.5, '  GPA' ) 

400  FORMATC '  ' ) 

RETURN 

END 
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SUBROUTINE  DUKP2(7AUGGR , TAUCGZ ,  MIACR , DRAGZ ,  QfDR ,  fkDZ , DELR , DELZ ) 
REAL  TAUCGR(0:6,0;6l),TAUCGZ(0:6,0:6l),DRAGR(0:6,C:6l), 

&  DRAGZ( 0:6,0:61 } ,D£LR,CELZ,L AMR, LAMZ,RLOC,ZLOC 
INTEGER  I,J,ENDR,EKDZ 
WRITE(7, 100) 

DO  10  I=0,EMDR 
RLOC=I»DELR 
DO  20  J=0,ENDZ 
ZLOC=J*DELZ 

LAMR=TAUGGR(I, J)/DRAGR(I,J) 

LAMZ=TAUGGZ( I , J)/DRAGZ( I, J) 

WRITE(7,200)  RLOC,ZLOC,TAUGGR( I , J) ,TAUGGZ(I , J) ,DRAGR( I , J) , 

&  DRACZ ( I , J ) , L AMR , L AMZ 

20  CONTINUE 

NRITE(7,30O) 

10  CONTINUE 

100  FORMAT ( '  1  MX,  'RLOC'  ,3X, 'ZLOC' ,3X,  'TAUGGR'  ,3X,  'TAUCGZ*  ,3X, 

&  ' DRAGR 4X, ' DRAGZ', 4X,' LAM  R\4X,’LAM  Z'/3X, 'CM' ,5* , 'CM' , 

&  nx^CDYN/cc'.sxj/eac-*)) 

200  FORMAT(2(1X,F6.3),6(1X,G8.2)) 

300  FORMAT ( '  ' ) 

RETURN 

END 


LIST  OF  SYMBOLS 


Symbol 


a 

b 

b-, 


P 

c  chem 


c  chem 
P 


Eign 


P9 


pg 


L 

n 

Pr 


Definition 


Units 


speed  of  sound  in  gas  cm/s 

burning  rate  coefficient  (cm/s )/( dyn/cm2 )n 

non-ideal  gas  equation  of  state  co-volume 
gas  specific  heat  at  constant  volume 
particle  specific  heat  at  constant  volume 


cm3/g 


erg/g*K 

erg/g*K 


gas-particle  interaction  drag, 
radial  component 


gas-particle  interaction  drag, 
axial  component 


gas  total  energy 
particle  total  energy 


gas  chemical  energy  released 
during  combustion 


particle  chemical  energy  released 
during  combustion 


particle  ignition  energy 

interphase  friction  factor 

convective  heat  transfer  coefficient 

gas  thermal  conductivity 

bulk  modulus  of  solid 

total  length  of  bed 

burning  rate  index 

gas  phase  pressure 

particle  phase  pressure 
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dyn/cm" 


dyn/cm3 

erg/g 

erg/g 


erg/g 


erg/g 

erg/g 


erg/s 'em2 *K 
erg/s 'cm'K 


g/cm*s2 


cm 


dyn/cm2 


dyn/cm* 


LIST  OF  SYMBOLS  (continued) 


Symbol 

Definition 

Units 

?% 

initial  gas  phase  pressure 

dyn/cm2 

\ 

initial  particle  phase  pressure 

dyn/cin2 

P  r 

gas  Prandtl  number 

0 

interphase  convective  heat  transfer 

erg/cm'**  s 

r 

radial  coordinate 

rp 

particle  radius 

ym 

"Po 

initial  particle  radius 

ym 

• 

r 

surface  burning  rate 

cm/s 

R 

total  bed  radius 

cm 

R* 

gas  constant 

erg/K*gmol 

Re 

Reynolds  number 

t 

time 

s 

T9 

gas  temperature 

K 

TP 

particle  temperature 

K 

\ 

initial  gas  temperature 

K 

TPo 

initial  particle  temperature 

K 

U9 

radial  component  of  gas  velocity 

cm/s 

UP 

radial  component  of  particle  velocity 

cm/ii 

W9 

axial  component  of  gas  velocity 

cm/s 

WP 

axial  component  of  particle  velocity 

cm/s 

z 

axial  coordinate 

h 

location  of  opening  in  outer  wall 

cm 

rc 

gas  production  rate  during  combustion 

g/cm^*s 
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LIST  OF  SYHBOLS  (concluded) 


Symbol  s 

<5 

A  r 
A  z 
Azc 
e 

m9 

u9o 

p9 

% 

P?o 

Pi 

n 

T 

i 


Def  Jn’t.ion 

finite  differencing  operator 
radial  distance  between  nodes 
axial  distance  between  nodes 
length  of  opening  in  outer  wall 
time  filter  coefficient 
gas  viscosity 
initial  gas  viscosity 
gas  phase  density 
particle  phase  density 
initial  gas  phase  density 
initial  particle  phase  density 
bulk  gas  density 
bulk  particle  density 
shear  stress  components 
porosity 


Units 


cm 

cm 

poise 

poise 

9/cm3 

g/cm3 

g/cm3 

g/cm3 

g/cm3 

g/cm3 

dyn/cm^ 


A 


« 


a 
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